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Abstract. This tutorial is a theoretical work, in which we study the physics of 
ultra-cold dipolar bosonic gases in optical lattices. Such gases consist of bosonic 
atoms or molecules that interact via dipolar forces, and that are cooled below 
the quantum degeneracy temperature, typically in the nK range. When such a 
degenerate quantum gas is loaded into an optical lattice produced by standing 
waves of laser light, new kinds of physical phenomena occur. These systems 
realize then extended Hubbard-type models, and can be brought to a strongly 
correlated regime. The physical properties of such gases, dominated by the 
long-range, anisotropic dipole-dipole interactions, are discussed using the mean- 
field approximations, and exact Quantum Monte Carlo techniques (the Worm 
algorithm). 
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1. Introduction 



In 1989, M. Fisher et. al. predicted that the homogeneous Bose-Hubbard model (BH) 
exhibits the Superfluid-Mott insulator (SF-MI) quantum phase transition T. In 2002 
the transition between these two phases were observed experimentally for the first time 
with cold atomic gases in the group of I. Bloch, T. Esslinger and T. Hansch [2]. The 
experimental realization of a dipolar Bose-Einstein condensate (BEC) of Chromium 
by the group of T. Pfau [3l HI [5], and the recent progresses in trapping and coohng 
of dipolar molecules by the groups of D. Jin and J. Ye [H [3 [8], have opened the 
path towards ultra-cold quantum gases with dominant dipole interactions. A natural 
evolution, and present challenge, on the experimental side is then to load dipolar BECs 
into optical lattices and study strongly correlated ultracold dipolar lattice gases. 

Studies of BH models with interactions extended to nearest neighbors had pointed 
out that novel quantum phases, like supersolid (SS) and checker board phases (CB) 
are expected [9l [TOl [TTl [12]. Due to the long-range character of the dipole-dipole 
interaction, which decays as the inverse cubic power of the distance, it is necessary 
to include more than one nearest neighbor to have a faithful quantitative description 
of dipolar systems. In fact, longer-range interactions tend to allow for and stabilize 
more novel phases. 

In this work we first study BH models with dipolar interactions, going beyond 
the ground state search. We consider a two-dimensional (2D) lattice where the dipoles 
are polarized perpendicularly to the 2D plane, resulting in an isotropic repulsive 
interaction. We use the mean-field approximations and a Gutzwiller Ansatz which 
are quite accurate and suitable to describe this system. We find that dipolar bosonic 
gas in 2D lattices exhibits a multitude of insulating metastable states, often competing 
with the ground state, similarly to a disordered system. We study in detail the fate of 
these metastable states, in particular what is their lifetime due to tunneling. Moreover, 
we find that the ground state is characterized by insulating checkerboard-like states 
with fractional filling factors v (average number of particles per site) that depend on 
the cut-off used for the interaction range. We confirm this prediction by studying 
the same system with Quantum Monte Carlo methods (the Worm algorithm). In this 
case no cut-off for the dipolar interaction is used, and we find evidence for a Devil' 
s staircase in the ground state, i.e. insulating phases which appear at all rational 
V of the underlying lattice. We also find regions of parameters where the ground 
state is a supersolid, obtained by doping the solids either with particles or vacancies. 
Recently [13], a complete devil' s staircase has been predicted in the phase diagram 
of a one-dimensional dipolar Bose gas. 

In this work, we also investigate how the previous scenario changes by considering 
a multi-layer structure. We focus on the simplest situation composed of two 2D 
layers in which the dipoles are polarized perpendicularly to the planes; the dipolar 
interaction is then repulsive for particles laying on the same plane, while it is attractive 
for particles at the same lattice site on different layers. Instead we consider inter- 
layer tunneling to be suppressed, which makes the system analogous to a bosonic 
mixture in a 2D lattice. Our calculations show that particles pair into composites, 
and demonstrate the existence of the novel Pair Super Solid (PSS) quantum phase. 

Moreover, we study a 2D lattice where the dipoles are free to point in 
both directions perpendicularly to the plane, which results in a nearest neighbor 
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repulsive (attractive) interaction for aligned (anti-aligned) dipoles. We find regions 
of parameters where the ground state of the system exhibits insulating phases with 
ferromagnetic or anti-ferromagnetic ordering, as well as with rational values of the 
average magnetization. Evidence for the existence of a Counterflow Super Solid (CSS) 
quantum phase is also presented. 

Our predictions have direct experimental consequences, and we hope that they 
will be soon checked in experiments with ultracold dipolar atomic and molecular gases. 

Although this paper reports on many novel results and predictions, it is written in 
a tutorial form. The reader will have a possibility to learn first the basic introduction 
to dilute Bose gases, and in particular dipolar Bose gases. We will explain very 
carefully various mean field methods that can be applied to describe approximately 
diluted dipolar Bose gases in 2D and 3D. These methods are elementary, but become 
more and more complex for the extended Hubbard models, and in particular for the 
models with infinite range interactions. In the chapters about multiple layers and 
up-down mixtures we will present and explain in detail methods of deriving effective 
low-energy Hamiltonians using second order degenerate perturbation theory. Last, 
but not least, we will present a very pedagogical introduction to QMC methods: path 
integral Monte Carlo method and the Worm algorithm. 

2. Dipolar Bose gas in optical lattices 

2.1. Optical lattices 

An optical lattice is an artificial crystal of light, resulting from the interference pattern 
of two or more counter propagating laser beams [IJ . The wavelengths Ai of the laser 
beams determine the spatial periodicity of the crystal; for example, two lasers of equal 
wavelengths A^; propagating along x but in opposite directions, produce a standing 
wave with an intensity pattern I{x) which is spatially periodic with periodicity Xx/2. 
An optical lattice can trap neutral atoms by exploiting the energy shift induced by 
the radiation on the atomic internal energy levels. 

The electric field E(r,i) = 2£'o cos (k • r — oj^t) of a monochromatic laser 
oscillating with frequency wl, interacts with a neutral atom, of spatial dimensions 
much smaller compared to the wavelengths of the light A^ = 2Tr/ki, {i = x,y,z), 
through the Hamiltonian 



where d = — e J2i is the electric dipole moment of the atom, the positions of the 
atomic electrons of charge e. With Hamiltonian (|2.ip . one can easily calculate the 
energy correction to the ground state of the atom, by means of perturbation theory. 
The first order correction vanishes because the dipole operator is odd with respect to 
space inversion (r^ — — r^), therefore the first non zero contribution is given by the 
second order correction 



i?int(t) 



d-E(r,<), 



(2.1) 



AE{r) 



-a(^i)(E(r,0')t 



(2.2) 



where 
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is the atomic polarizability |15[ 116] . and {■■■)t denotes a time average over one 
oscillation period of the electric field [J. In the last expression the energies in the 
denominators are the unperturbed energies of the atom, where g is the ground state, 
the sum runs over all excited states 7, and e is the unit vector in the direction of 
the electric field. In a typical experiment the laser light is far off resonance, which 
means that the laser frequency is close to one of the unperturbed excited states (e.g 
Ee = huje), but does not induce any real transition. In such a situation, one can take 
only the smallest of the denominators (|2.3p . and the polarizability becomes inversely 
proportional to the laser detuning from resonance fiA = Hujl — {E^ ~ Eg) 

|(e|d.g|g)|^ 

^^A ■ (^-^^ 

In this situation, the energy shift is given by 

AE{r) ^ -iaK)(E(r,t)2), « (2.5) 

where /(r) is the intensity of the laser. In the dressed atom picture, the energy shift 
(|2.5|) is interpreted as an effective potential Vopt{j^) = Ai?(r), that follows the spatial 
pattern of the laser field intensity, in which the atom moves. In this picture, the atom 
then feels a force 

Fdipoio = -VFopt(r), (2.6) 

that attracts it towards the regions of high intensity for the so called red-detuned 
lasers (i.e. A < 0), while a blue-detuned light (i.e. A > 0) pushes the atom out of 
the regions of high intensity. In the literature this force is called the dipole force, as 
it is the resulting interaction of the induced atomic dipole moment with the spatially 
varying electric field of the light. Note that in order to reduce heating caused by 
inelastic scattering, i.e. photon absorption and spontaneous emission processes, a 
large detuning is required because the photon scattering rate scales as /(r)/A^. In 
the limit of large detuning an optical lattice is therefore non-dissipative, which makes 
it a basic tool to manipulate cold neutral atoms. 

For example, the simplest case of a one-dimensional lattice is obtained by the 
superposition of two lasers propagating in opposite directions, with electric fields 
linearly polarized, say in the z direction, and given by 

Ez{x, t) = 2Eq cos (kxX — ujLt) + 2Eo cos {—k^x — LOi,t) 

= 4i?o cos (fca;^) cos(ajLt). (2.7) 

The time average over one period of oscillation of the electric field gives then 
{Ez{x,t)'^)t — 2EoCOs'^{kxx), which yields to the spatially varying optical potential 

Vopt{x) = Vb,^ cos^ik^x), (2.8) 

with periodicity Xx/2 = n/kx, and Vq^x = 2Eoa{ujL) from Eq. p.4p . The 
generalization to the two dimensional (2D) or three dimensional case (3D) is 
straightforward (see e.g. ^15]). For example in Fig. [T] two different geometries are 
shown. 

X more specifically (■■■}* = 7 /q ■■■ dt where t = nn/uji^, n = 1, 2, ■ • • 
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Figure 1. Picture of an optical potential, (a) 2D square lattice of quasi ID traps; 
(b) a 3D cubic lattice, picture courtesy of I. Bloch I14| . 



2.2. Theory of dilute Base gases 

In this section we recall some basic theory of a dilute gas of neutral bosonic particles 
at temperature T well below the degeneracy temperature. The type of particles we 
consider here can be atoms or molecules. 

For a dilute gas, the interparticle separation (typically of the order of 10^ nm for 
alkali atoms [15 J is an order of magnitude larger than the length scales associated 
with the atom-atom interaction. In other words, a dilute gas of density n is a very 
rarefied gas in which the "spatial extension" of an atom is much smaller than the 
average volume per particle n~^. Because of this condition, the two-body interaction 
dominates the physics while three-body or more are very unlikely and essentially not 
important. The two-body interatomic potential ^(r) depends on the type of particles 
one considers, the relative distance between the atoms r = ri — r2 and their internal 
states. For alkali atoms, the potential is strongly repulsive for small atomic separations 
while for large atomic distances it is dominated by the van der Waals attractions that 
decay as — Cg/r^, where the coefficient Cg depends on the atomic species. 

Here we will consider only elastic scattering, where the internal states of the 
two atoms do not change in the collision process. If the temperature of the gas is 
very low, i.e. T — > 0, the kinetic energy of the particles is very small compared to 
the centrifugal barrier and only s-wave scattering takes place. Therefore, the only 
important parameter is the scattering length given by 



with m being the mass of the atoms. This quantity has the dimensions of a length, 
and for Os > has the physical interpretation of the radius the atoms would have if 




(2.9) 
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they were considered to be perfect billiard balls. The condition for the diluteness of 
the gas then reads 

n|as|3<l (2.10) 

where n is the density of the gas and njagp is called the gas parameter. One can 
invert expression (j2.9p and think of an effective contact interaction between the two 
particles proportional to the scattering length, and given by 

yeff(r)=ff<5(3)(r), (2.11) 

where g is defined as 

5 = , 2.12 

m 

and S is the Dirac delta function. Note also that since the effective interaction depends 
only on the scattering length, it is repulsive (attractive) for positive (negative) a, and 
it can be dynamically modified for example in alkali atoms just by varying an external 
magnetic field near a Feshbach resonance. 



2.2.1. The Gross- Pitaevskii equation — The quantum state of a gas of N particles is 
described by the many-body wavefunction 'J'(ri,r2, ...jTn), and the time evolution of 
the system is determined by the Schrodinger equation. In a BEC, one can describe the 
dynamics of the condensate just through the Gross-Pitaevskii (GP) equation [15] 
given by 

z;i^vE'o(r, t) = + ^cxt(r) + 5l^o(r, t)\^^ ^oir, t), (2.13) 

where ^'o(r,i) is the BEC wavefunction, also called the order parameter. 
The interaction between particles has been taken into account in a mean-field 
approximation by the term g\'^o{r,t)\'^ , where g is given in Eq. 12.121 Vc^t{^) is 
an external trapping potential, and the order parameter is normalized to the total 
number of particles, i.e. N = J d"^r|^'o(r, t)p. Equation ()2.13|) was independently 
derived by Gross and Pitaevskii in 1961, it is one of the main theoretical tools 
for investigating dilute weakly interacting Bose gases at low temperatures, and it 
has the typical form of a mean field equation where the order parameter must be 
calculated in a self-consistent way. The GP equation has proven to be a very useful 
tool to describe the physics of weakly interacting Bose-Einstein atomic condensates 
in the early ages of this field. With this formalisms, and its extension to include 
small fluctuations given by Bogoliubov theory, one can describe accurately, among 
others, the collective excitations of the systems, the response to rotations including the 
formation of vortices, the propagation of sound, the presence of dynamical instabilities. 
Generally speaking, the GP treatment is well suited in the regime of full coherence, 
when a single macroscopically occupied matterwave correctly describes the system. 
At the end of the '90, few years after the creation of the first alkali BECs in the 
lab, the need of "going beyond GP" started to be very strongly felt, due to the 
theoretical interest and experimental possibility of going into the strongly correlated 
regime. In fact, the presence of strong interactions, strong rotations and/or special 
trapping potentials can limit the validity of the GP equation. For instance a strong 
confinement in one or two dimensions can reduce the system to an effectively 2D or ID 
one. A strong rotation combined with interactions can lead to quantum Hall physics. 
Also the presence of a deep optical lattices, when the combined effect of interactions 
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and trapping potential leads to a "fragmentation" of the condensate, requires more 
sophisticated descriptions. 

In this tutorial we are interested in describing the physics of Bosons trapped in a 
periodic optical potential (Vopt) and eventually also confined in a magnetic harmonic 
trap (Vho), the total external field being given by the sum 

K>xt(r) = Kpt(r) + 14o(r) 

= Vo,^cos^{hri) + ^m ^ cofrl (2.14) 

i—x,y,z i—x^y^z 

where {Vq.x, Vb,y, Vb,z) is the depth of the optical lattice in the three spatial directions 
and (wx, Wy , Wz) the frequencies of the harmonic trap. In order to describe the physics 
of Bosons trapped in the potential (I2.14[) . we need to "go beyond" the GP equation, 
and we will devote the following sections to this purpose. 



2.2.2. Bose-Hubbard model — The starting point of our discussion is Hamiltonian 
(j2.15l) . written in the second quantization formalism in terms of the creation and 
annihilation operators for Bosons, ip^ir) and '!/'(r) respectively, and given by the 
expression 



2m ' ' 2 



^(r), (2.15) 



where the first term in square brackets is the kinetic energy, Vi^xt (r) — Kpt (r) + V^o (r) 
is the external trapping potential (j2.14p and we have used the simplified contact 
interaction (|2.1ip . We work in the grand canonical ensemble and the chemical 
potential fx fixes the total number of particles. Additionally, we assume the harmonic 
confinement to change on a scale larger than the one of the optical lattice, such that 
we can consider the effect of the magnetic trapping to be constant over a single site 
of the lattice. 

In this formalism, the field operators can be written in the basis of single-particle 
wave functions {$„(r)}„, where ti is a complete set of single particle quantum numbers 

(2.16) 

^t(r)=^co;(r)at, 

n 

with a|j and a„ being the creation and annihilation operators for the mode n, i.e. 
aj'jn) = y/n + l\n + 1) and d„|n) = ^/n\n — 1). Also, the field operators satisfy the 
usual commutation relations for Bosons 

oo 

[V'(r), ^t(r')] = ^ $„(r)$:(r') = ^^(r ~ r'), ^ ^ 

^0 (2.17) 

[7^(r),^(r')] = [7^^(r),^t(r')] ^0. 

It is well known |18| , that the spectrum of a single particle in a periodic potential 
is characterized by bands of allowed energies and energy gaps, and the single particle 
wave functions are described by Bloch functions $Qk(r) with band index a and quasi- 
momentum /ik. Alternatively, there exists a complementary single-particle basis given 
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by the Wannier functions [IHIIII] ^^(r — R^), where Rj is a lattice vector pointing at 
site i and (r) are defined as the Fourier transform of Bloch functions 

«;„(r) = ^^e-k.r<j,^k(r), (2.18) 



k 

where Ns, is the total number of sites in the lattice. The Wannier functions form a 
complete orthonormal set, so one may write the field operators (|2.16p as 

(2.19) 

^t(r) = *«k(r)aL = E " 

Wannier functions are useful in the case of deep optical lattices where tight binding 
approximation apply. The big advantage of using Wannier functions ^^(r ~ Ri) is 
that they are localized and centered around the lattice site pointed by R.;. 

If the temperature of the system is low enough, and the interactions between the 
particles is not sufficient to induce transitions between the bands, one may restrict only 
to the first Bloch band because the particles have insufficient energy to overcome the 
gap that separates the first band from the others. This amounts to keep in (|2.19p only 
the lowest of the a indices, which we omit for simplicity of notation in the following. 
Therefore the Hamiltonian (|2.15p becomes 

H = - J^j a\aj + E ^-|^ a\a\akai - ^ a\aj, (2.20) 
where the quantities in the sums are given by 



Jij = - / d'^rw*{r - Ki 



2m 



+ K,pt(r) 



w(r - Rj) (2.21) 



U^jki=9 J dW(r-R,)w*(r-R,Xr-Rfe)w;(r-RO (2.22) 

^i^J = y"dW(r-R,)[/i-Ko(r)]w(r-R,). (2.23) 

The Wannier functions are localized on the lattice sites, the deeper the lattice the more 
localized they are. For a sufficiently deep optical potential, in Eq. (|2.22p and (|2.23p 
the dominant contributions are given by Una and /in. For the kinetic part (|2.2ip . 
there is a constant contribution given by Ja and due to the presence of the derivative 
in the integration, there is also a positive matrix element for nearest neighboring 
sites Jij > 0. The two situations are qualitatively shown in Fig. ^ where we have 
approximated the Wannier functions with two Gaussians respectively localized at site 
i and j of the lattice. However, we stress that the picture provided by Gaussian 
functions is only qualitative. In fact, in order to be quantitatively correct, one needs 
to calculate the proper matrix elements with Wannier functions. 

With the above considerations, we can now write the celebrated Bose-Hubbard 
Hamiltonian in the form 

Hbh = E + ^ E ~ -'^) ~ E f^''^'' (2.24) 

(ij) i i 
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Figure 2. (a) Two Gaussians localized on neighboring sites i and j of an 
optical lattice having negligible overlap, (b) The first derivatives of the Gaussian 
functions instead, show a negative overlap in the region indicated by the arrow, 
which leads to a positive matrix element Jij > 0. 

where (ij) indicates sum over nearest neighbors, the tunnehng coefficient J — Jij — Jji 
for hermiticity, the on-site interaction U = g f d'^r \w{r)\'^ , fii = a\ai is the number 
operator at site i, and we have neglected J a since it gives a constant contribution for 
each site. The harmonic confinement, since it is assumed to be constant across one 
lattice site, has been taken into account in the chemical potential as 

M» =M- ^"1^2 . (R, -Ro)2, (2.25) 

where Ro is the center of the harmonic trap with frequencies given by w = (wx, Wy , Wz) 
in the three directions. The second term on the right hand side of Eq. (|2.25l) is 
practically a chemical potential that differs from site to site and it is often called the 
local chemical potential. 

For a one dimensional optical lattice Vo-pt{x) = Vbsin^(fca;) with wavevector 
k = 2it/\, Fig. [3] shows both the on-site interaction U (solid line) and the tunneling 
coefficient J (dashed line) as a function of the optical lattice depth Vq, where all 
the quantities are measured in terms of the recoil energy Efi = S^fc^/2m, that is the 
energy acquired by the atom after absorbing a photon with momentum hk. The lattice 
parameters U and J were calculated numerically in e.g. |20| for different values of Vq. 
From Fig. |3] (b) , it is clear that it is possible to change the tunneling coefficient J over 
a wide range, going from a situation of practically isolated lattice sites at Vq — 25Eii 
up to a regime in which particles can tunnel from site to site at Vq = 5Efi, only by 
changing the optical potential depth by a few tens of recoil energies, and leaving the 
on-site interaction U practically unchanged. 

2.3. Dipolar Base gas 

2.3.1. Properties of the dipole-dipole interaction — Two particles 1 and 2 in a three 
dimensional space, at relative distance r and with dipole moments along the unit 
vectors ei and 62 as in Fig. |4](a), interact through the dipole-dipole interaction such 
that their interaction energy is given by 

An r° 

where r = |r|, and C/dd(r) = C/dd(— r). The dipolar coupling constant Cdd is different 
for particles having a permanent magnetic dipole moment fj,, and for particles having 
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Figure 3. (a) Schematic representation of a ID optical lattice, where s is a 
local chemical potential; (b) scaled on-site U (solid line) and tunneling coefficient 
J (dashed line) dependence on the optical potential depth Vg. The on-site 
interaction is multiplied by a/as{^ 1), where a = A/2 is the lattice period and 
as is the s-wave scattering length for atoms of equal mass m. Figure courtesy of 
D. Jaksch [20]. 



a permanent electric dipole moment d, and is respectively given by 

_ / MoAi^ magnetic , , 

^'^^"t dVeo electric, 

where fiQ is the vacuum permeability, and Eq is the vacuum permittivity. 

The dipole-dipole interaction (j2.26p has a long-range character; this is because at 
large distances it decays as Udd ^ ^/r^, contrary to the typical van der Waals potential 
that behaves like U^dw Also, from (j2.26p it is easy to see the anisotropic 

property of this interaction; for polarized atoms, i.e. all dipoles pointing in the same 
direction, the interaction reduces to 

C/d.(r) = ^i:i^, (2.28) 

where is the angle between the dipole and the relative distance of the particles, as 
in Fig m (b). The interaction is repulsive for 9 — 7r/2 as the example of Fig 2] (c), 
and attractive for 9 = as shown in Fig |4] (d) . The situation is reversed for anti- 
parallel dipoles, where a minus sign appears in front of Eq. (j2.28p . and therefore the 
interaction is attractive for 9 = tt/2 while 9 = gives rise to repulsion. 

The scattering properties of ultracold atoms, in the simple case of isotropic van 
der Waals interactions, are entirely described by the s-wave scattering length and the 
potential can be replaced by the effective contact interaction (j2.1ip . In the presence 
of a dipolar interaction as (I2.26p . because of its long range (decay as 1/r'^) and 
anisotropic character (strong dependence on the relative angles between the dipoles), 
all partial waves contribute to the scattering problem and also partial waves with 
different angular momenta couple with each other. While for Fermions, replacing 
the real potential (|2.26p with an effective dipolar interaction as (|2.1ip is reasonable 
[21], for Bosons this is not obvious, and in recent years it has been the subject of 
intensive studies [HI [2H Hi]. In the presence of an optical lattice, it has been 
recently argued [37] that in a ID geometry, replacing the real dipolar potential with 
an effective interaction as (|2.1ip is reasonable as long as the optical lattice is shallow 
enough. However, in the most general case it is necessary to account for the full 
expression of the dipole-dipole interaction potential (|2.26p . 
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Figure 4. (a) Two dipoles, 1 and 2, directed along unit vectors ei and e2 and 
separated by a distance r. (b) Polarized dipoles, for which the interaction depends 
on the angle 8 between the direction of the dipoles and the interparticlo separation 
r. This results in a repulsive interaction for 8 = tt/2 (c), and attractive for 6 = 0,it 
(d). 



2.3.2. Polarized dipoles in anisotropic harmonic traps — We now move to the 
description of a BEC of polarized dipoles, pointing along the z axis. For polarized 
dipolar BECs, due to the anisotropy of the dipolar interactions, the geometry of 
the trapping potential plays a fundamental role, first in determining the spatial 
distribution of the density, and second in the stability of the gas. 

Qualitatively, there are two extreme scenarios depending on the shape of the 
confining potential, shown in Fig. [S] (i) for a cigar-shaped trap elongated along the z 
axis, i.e. with an aspect ratio between the axial lJz and radial frequencies ujp — uj^ — Wj, 
given by A = ujz/ujp <^ 1, the density is mainly distributed along the polarization axis 
and the effect of dipole-dipole interaction is mostly attractive, which might lead to an 
instability of the gas even in the presence of a weak repulsive contact interaction; (ii), 
for a pancake-shaped trap, which is strongly confining along the z axis, i.e. A 3> 1, 
the dipolar interaction is mostly repulsive and the BEC is always stable for repulsive 
contact interactions and might be stable even for attractive contact interactions. In 
an intermediate situation in which the confining potential is perfectly spherical, the 
density distribution is then isotropic and the dipole-dipole interaction averages out to 
zero, which leads to a stable BEC for repulsive contact interactions. One can switch 
between one or the other scenario, just by adjusting the frequency of the confining 
potential along the z axis with respect to the axial x and y, and therefore it is natural 
to expect that for any given A there is a critical value for the scattering length Ccrit 
below which the BEC is unstable [28] . 

One can quantitatively describe the above scenarios starting from the Hamiltonian of 
the system, which in the presence of the dipole-dipole interaction (|2.28p reads 

H= j A\i;\v) _:^ + 1/,t(r) + |7A^(r)7A(r)-A* ^(r) 

+ i / d\id\2i^\v^)i;\v2)UAa{vi - r2)Vi(ri)^(r2). (2.29) 
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(a) 



(b) 



Figure 5. Polarized dipoles in anisotropic harmonic potentials, (a) in a cigar 
shaped trap elongated in the direction of polarization, the resulting dipolar 
interaction is attractive, and (b) in a pancake trap with a strong confinement 
in the direction of polarization, the dipolar interactions are repulsive. 



With the same approximations used to derive the Gross-Pitaevskii equation, one can 
write the Boson field operator ?/;(r) = 5'o(r) + (5?/;(r) as a sum of a classical field \I/o(r), 
the condensate wave function, plus the non condensate component dip(j) [17]. By 
neglecting the fluctuations Sil;{r), one can calculate the energy of the BEG given by 



2m 
^l*o(r) 



iVvfoWI^ + KxtWI^oWl^ 



{/dd(r-r')|*o(r')rdV 



2 

dV, 



*o(r)r 



(2.30) 



where the macroscopic wavefunction is normalized to the total number of particles 
N. Within a variational Ansatz, we assume the condensate wave function to be a 
Gaussian of axial width and radial width ax = Uy = cTp, normalized to the total 
number of particles N, namely 



*o(^;,/o) = 



N 



exp 



1 



Ho 



P 



(2.31) 



where aho — \/h/ (muj) is the harmonic oscillator length with average trap frequency 
uj = {LUpCOzY'^. Therefore, inserting Ansatz (|2.3ip into the energy functional Eq. 
(|2.30p , after integration we find the energy of the BEC to be a function of the widths 
of the Gaussians, namely 



with the kinetic energy 

Nhu) ( 2 

Eun — -. + 

the potential energy due to the trap 



trap 



Er. 



E, 



dd, 



EtYap — 



Nhui 



2<T +X <Jz' 



4A2/3 

the contact interaction energy given by 

hu) 1 



E, 



contact ^ 



/27raho f^pf^z 

and the contribution coming from the dipolar term 

toodd 1 



Edd — — 



27raho (^pCTz 



(2.32) 
(2.33) 

(2.34) 
(2.35) 

(2.36) 
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We have introduced the dipolar length a^^ = -§7^^, with C^d given in Eq. (|2.27p . 
which measures the absolute strength of the dipolar interaction, k = (Jp/a-^ is the 
aspect ratio of the density distribution, and the function / is given by 



1 + 3K^artanh-\/l — 

^('^^^T^^ (1 - ^2)3/2 • (2.37) 

While the integrals needed to obtain (|2.33I2.34I2.35P are easy to calculate since they 
contain only Gaussian functions and their derivatives, the integral to get (|2.36p is not 
straightforward due to the presence of the dipolar potential [/dd(ri — r2). See section 
l2.3.3l for more details. In the left panel of Fig. [Bl we show the behavior of the function 
/(k) as K is continuously varied from k — 10^^ to k = 10^. The function takes the 
asymptotic values of /(O) = 1, /(oo) — ~2, and it vanishes for k = 1, which implies 
that for a spherical density distribution the dipole-dipole mean-field interaction (j2.36p 
averages out to zero. Therefore we notice that it is possible to control the strength 
and the sign of the mean-field dipolar interaction just by adjusting the aspect ratio A 
between the axial and the radial frequencies of the confining trap. The total interaction 
energy is provided by the sum of the contact (j2.35p plus the dipolar interaction energy 
((236)) . given by 

E.n. = ^^(^-m). (2.38) 



27raho o'pO'z \add 

The stability of the gas requires a repulsive interaction E-mt > 0, which leads to the 
condition 

— - f{n) > 0, (2.39) 

Odd 

and can be adjusted ad-hoc by changing the frequencies of the trap in the three 
directions. 




0.01 0.1 1 10 100 1" "1 1 1" "1 

K Trap aspect ratio A = uJz/ujp 



Figure 6. (Left panel), ft dependence of the /(k) function that appears in 
the mean-field dipolar interaction. (Right panel) Stability diagram of a dipolar 
condensate: the thin line is the solution for a^^n^{X)/ao calculated with the 
Gaussian Ansatz 1 12. 311 1. where ag is the s-wave scattering length, while the thick 
line is the numerical solution of the GP equation 1291 . The dots with error bars 
are experimental data [31) . 
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To determine the stability threshold acrit(A), one needs to minimize the energy (|2.32p 
with respect to the variational parameters ap and for fixed values of N, A and w. 
The results are summarized in the right panel of Fig. IH] as a thin line, while the thick 
line represents more accurate results calculated from solving numerically the Gross- 
Pitaevskii equation |29j . The dots with error bars correspond to experimental data 

m- 



2.3.3. Mean- field dipolar interaction in a spherical trap — In order to calculate the 
mean-field dipolar interaction energy (|2.36p , we insert the Gaussian Ansatz (|2.3ip into 
the the second of the integrals (|2.30p . and we get to the expression 



Edd^^^J dVidV2p(ri)p(r2)C/dd(ri -ra), 



(2.40) 



with p(r) = |^'o(r)p being the condensate density at r. The last integral can be 
simplified by means of the convolution theorem j28[ [50] which states 



J dV2;7dd(ri - r2)p(r2) = {c/dd(k) p(k)} 



(2.41) 



where t/dd(k) and p(k) are the Fourier transform respectively of the dipole-dipole 
potential and the density. J-~^ indicates the inverse Fourier transform, and using its 
definition we can write 



dd 



i / dVip(ri) 



1 



1 



(27r)3 
d3fcC^(k) p\k), 



d3A:C/dd(k)p(k)e 



zk-ri 



2(27r)3 

where in the last step we have used the relation p(k) = p(— k). 

The Fourier transform of the dipole-dipole interaction (I2.28P is given by 



f/dd(k) = / d3r[/dd(r)e 



Cdd(cos^7 - 1/3), 



(2.42) 



(2.43) 



where 7 is the angle between k and the polarization direction, and Cdd is given by 
expression (|2.27|) 28 . In order to evaluate the integral of Eq. (|2.42p we need to insert 
the condensate wave function, which in the simple case of isotropic potential (cr^ = CTp) 
becomes a product of three Gaussian distributions with equal widths a. Therefore the 
Fourier transform of the condensate density is readily calculated as 



N 



(VTTCraho)^ 



d^re-'^-^e 



exp 



2„2 



a a 



ho 1,2 



(2.44) 



This expression has to be inserted into the integral (|2.42l) . which can be easily 
evaluated in polar (r, 7, 1^9) coordinates giving 

NCdd 



Edd — 



■J sin7d7d(/3fc^dA:(cos^ 7 — 1/3) exp 



^hoi 



= 7VCdd27r / dfcfc^exp 
0, 



^2„2 



ho 



2 

1/3) 



(2.45) 

where we have performed the change of variable x = cos 7. The generalization to 
anisotropic density distributions is mathematically more demanding but in principle 
straightforward, and leads to Eq. (I2.36p . 



§ remember, 7 is the angle between the polarization axis z and k. 
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2.3.4- Extended Bose-Hubhard model — As in section [2.2.21 we expand the field 
operators in the basis of Wannier functions (|2.19p , and we keep only the lowest index 
corresponding to the first Bloch band. Within this approximation the first line of Eq. 
(|2.29p leads to the Bose- Hubbard Hamiltonian (|2.24p . Instead the dipolar term gives 
rise to a further contribution 

Hdd^^J dVidV2V3^(ri)^^(r2)C/dd(ri - r2)^(ri)^(r2), (2.46) 

which after expansion of the field operators in the basis of Wannier functions, becomes 

Hdd^ (2.47) 

i,j,k,l 

and the matrix elements Vijki are given by the integral 

V.jki = J dVidV2U.*(ri - R0u'*(r2 - R,)[/dd(ri - r2)u;(ri - Rk)w{r2 - R/).(2.48) 

The Wannier functions are centered at the bottom of the optical lattice wells with 
a spatial localization that we assume to be a. For deep enough optical potentials 
we can assume a to be much smaller than the optical lattice spacing d, i.e. a <^ d. 
In this limit, each function w{r — R^) is significantly non-zero for r ^ R^, and the 
integral (|2.48|) is significantly non-zero for the indices i = k and j — I. Therefore there 
are two main contributions to the integral (|2.48p : the off-site matrix element Vijij 
corresponding to k = i^j = l, and the on-site Van when all the indices are equal. 
Below we will explain the physical meaning of these two contributions. 



Off-site — The dipolar potential C/dd(ri — V2) changes slowly on the scale of ct, 
therefore one may approximate it with the constant C/dd(R.j: — Rj) and take it out of 
the integration. Then the integral reduces to 



Vijij 



C,„(R,-R,, / d=n IMr, - R.)l7dV. - R,,r , (2.49) 



which leads to the off-site Hamiltonian 

2 



In the last expression Vij — C/dd(R-j: ~ R-i)i "-i — o^iO-i is the bosonic number operator 
at site i, and the sum runs over all different sites of the lattice. 



On-site — At the same lattice site j, where |ri— r2| ^ cr, the dipolar potential changes 
very rapidly and diverges for |ri — r2| — > 0. Therefore the above approximation is not 
valid any more and the integral 

V^u^ = j dVidV2p(ri){/dd(ri -r2)p(r2), (2.51) 

with yo(r) = |w(r)|^ being the single particle density, has to be calculated taking into 
account the atomic spatial distribution at the lattice site, similarly to what has been 
described in Sec. 12.3.21 j|. We have already encountered this kind of integral in Sec. 



11 Since Ri is a constant, we have renamed the variables as — = for u = 1,2. 
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12.3.21 and we have seen that, a part from a factor of 2, the solution can be found by 
Fourier transforming, i.e. 

Vuu = J d^kU^d(k.)p^{k). (2.52) 

Which leads to an on-site dipolar contribution to the Hamiltonian of the type 

i 

The extended Bose-Hubbard Hamiltonian is given by the sum of the Bose-Hubbard 
(|2.24p and the dipolar Hamiltonians calculated above, leading to the expression 

^cBH = "i^i + ^ ^ ^ 5Z ~2^^* "j' (2-54) 

(ij) i i i^j 

where U is now taken into account as an effective on-site interaction 

U = gJ d^r\wir)f + -^ J d^kU^^ik) p{k), (2.55) 

which contains the contribution of the contact potential, with g given in Eq. (|2.12p . 
plus the dipolar contribution coming from (I2.53p . Approximating each lattice site 
with a tiny harmonic trap, and approximating the atomic density distribution with 
Gaussians, U looks like Eq. p.38p . and one can see that the resulting on-site 
interaction can be increased or decreased by changing the lattice confinement. 



3. Hubbard models: theoretical methods 



3.1. Superfluid— Mott insulator quantum phase transition in the Bose-Hubbard model 
Consider the Bose-Hubbard Hamiltonian as derived in Sec. ()2.2.2p . 

with a uniform chemical potential /i, and a total number of Bosons given by the 
expectation value of the operator N = '^^fii. There are three parameters in this 
Hamiltonian, namely J, U and /i, but it is a convention to reduce the analysis of the 
phase diagram of Hbh to the ratio of two of them over the third one, e.g. J/U and 

t^/u. 

The ground state of Hamiltonian (j3.ip is easily understood for two opposite 
regimes of parameters: (i) for shallow lattices, i.e. U/J <^ 1, the system is in a 
gapless superfiuid phase (SF) characterized by on-site density fluctuations and the 
particles delocalized over the whole lattice; (ii) for deep lattices, i.e J/U <C 1, and 
commensurate filling, on-site density fluctuations are completely suppressed, each site 
is occupied by an integer number of atoms n, and the ground state is a product of 
single-site Fock states 

\GS) ^\n,n,---). (3.2) 

This filling is energetically favorable in the range of chemical potential n—l< ji/U < 
n. The system is gapped and incompressible, as beautifully explained in the famous 
paper of Fisher et. al. [1], and it is called a Mott insulator MI(n). For small values 
of J/U the MI(n) phase persists in a closed and finite area of the J vs. /i plane 
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[T], which is called the Mott lobe for MI(n). The larger J/U value of the lobe is 
called the tip of the lobe or also critical point {J/U)c- The critical point changes 
with the dimensionality and geometry of the system. In Fig. [7] we plot the first 
n = 0,1,2,3 insulating lobes, calculated for an infinite optical lattice within the 
mean-field approximation, which will be discussed in Sec. 13.2.21 The thick black lines 
enclose the lobes and mark the boundaries between the MI and SF phases. Outside 
the insulating lobes, the phase is SF. The colored lines of Fig. [7l[a) indicate a contour 
plot of constant fractional density, while the thick black lines departing from the tip 
of the lobes and extending into the SF region, correspond to an integer value n of the 
density. 





0.08 0.16 0.24 0.04 0.16 0.24 
zJAJ zJ/U 



Figure 7. Mean-field phase diagram of the Bose Hubbard Hamiltonian. (a) 
Contour plot of the density per site; (b) contour plot of the order parameter (see 
Sec. 13.211 . MI(n) indicate a Mott insulating phase with fixed n atoms per site. 



We will derive the mean- field Mott insulating lobes of Fig. [7] in a more rigorous 
way in Sec. I3.2.2[ but for the moment we just list the critical points {J/U)c, 
for the n = 1 lobe, that have been estimated with different methods and for 
different dimensions of the lattice. In one dimension, the critical point has been 
estimated to be {J/U)c — 0.29 |31] using Density Matrix Renormalization Group 
calculations (DMRG). In two dimensions, with quantum Monte Carlo calculations, 
the critical point has been estimated to be {J/U)c — 0.061 while in the 

three dimensional model the location of the critical point has been estimated with 
perturbative expansions [34], and quantum Monte Carlo simulations |35| to be at 
{J/U)c — 0.034. In the next section we will derive the mean-field lobes. 
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3.2. The Gutzwiller mean-field approach 

The Gutzwiller mean-field approach provides an approximated many-body wave 
function for Hubbard-type Hamiltonians, of the form 

I*) =n (3-3) 

i n—Q 

where \n)i represents the Fock state of n atoms occupying the site i, n^ax is a cut off 
in the maximum number of atoms per site, and is the probability amplitude of 
having the site i occupied by n atoms. The probability amplitudes are normalized to 
unity E„ = 1- The Gutzwiller wave function p.3p has been extensively used 

in the literature [201 [36l |37l |38] , it predicts that there exists a critical value {J/U)mf 
in a given range of /x, below which the ground state is a product of single Fock states 
fn^ = Sn,n with cxactly n particles per site, as p. 21) . Moreover, for .J/U > {J/U)mf the 
Gutwiller Ansatz predicts a superfluid ground state with fluctuating on-site particle 
number. The Gutzwiller critical point, for n = 1, is found to be (J/C/)mf = 1/5. 8z 
[55] . where z = J2{j)i^ is the number of nearest neighbor connections at each site of 
the lattice. In table [T] we show the comparison of the critical points predicted by the 
Gutzwiller Ansatz for different dimensions of the system, with the more precise ones 
discussed in Sec. p.l[) . From the comparison, one can deduce that the Guzwiller is 



Table 1. Comparison of the Gutzwiller critical points {■J/U)^{ with the more 
precise, up to now, critical points {J/U)c, for different dimensions D of the system. 



D 


z 


{J/uu 


{J/U)c 


1 


2 


0.0862 


0.29 


2 


4 


0.0431 


0.061 


3 


6 


0.0287 


0.034 



unsatisfactory for ID systems {z — 2) while it is satisfactory for a 3D one (z = 6). 
Also, in the limit of J/U — ?> oo the difference between the Gutzwiller predictions and 
the exact results are negligible [39 . Summarizing, the Gutzwiller predictions are exact 
in the two limiting cases oi J/U — J> 0, and J/U — > oo, while for intermediate cases 
the performance of the Gutzwiller approach strongly depends on the dimensionality 
of the lattice, since it does not correctly account for the quantum fluctuations at the 
phase transition. 

Two important quantities are the order parameter ipi = (^'|ai|^'), namely the 
expectation value of the Bosonic annihilation operator at the i-th site of the lattice, 
and the density fluctuations at the site i given by Sui = — (vI/|77,j|vl/^2_ gy 

using the Gutzwiller wavefunction p.3p one gets the following expressions for the order 
parameter 

(3.4) 

n 

and for the density fluctuations 

5n.=^n2|/«p-[^n|/Wp]'. (3.5) 
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The order parameter tpi together with the density fluctuations Sui describe the phase 
of the system at the site i of the lattice: in the Mott phase ipi = 0, and density 
fluctuations are suppressed Srii = 0, while the superfluid phase is characterized by 
(fi ^ and presence of density fluctuations Sui ^0. In the uniform system, the 
lattice is translationally invariant and therefore all sites are self-similar, which means 
that only one site is sufficient to determine the phase of the whole system. In Fig. [TJb) 
we plot the absolute value of the order parameter ifi for such a system, in the J/U vs. 
^/U plane. The colored lines outside the insulating lobes correspond to a contour plot 
of constant non-zero value of tpi typical of the SF phase. Instead, in a non-uniform 
system as it is in the presence of an external confining harmonic potential, different 
phases can coexist. As an example, in Fig. |8]we plot the density of the ground state 
(a), the order parameter at each site (b), and the density fiuctuations (c) of a 2D 
lattice in the presence of a confining harmonic potential. Notice that the MI phase 
at the center of the harmonic trap (0,0) is surrounded by a ring of SF phase, in a 
wedding-cake like structure, as first discussed in [3^. The ground state of Fig. [5] was 
obtained within the mean-field approximation through the imaginary-time evolution 
technique, that will be discussed in Sec. 13.2.11 



(a) W (c) 




10 -10 10 -10 10 -10 



Figure 8. Mean-field ground state of a 2D optical lattice, calculated for 
J = 0.025C/, and fii = 0A5U - fii^, where U = 8 X IQ-^U/h is the frequency 
of the harmonic oscillator confinement, (a) the vertical axis shows the value of 
the density at each site of the lattice, (b) the corresponding absolute value of the 
order parameter, and (c) on-site density fluctuations. 



In the presence of dipolar interactions, as we shall see later on, it is also necessary 
to account for non-uniform quantum phases, because even in the uniform system 
the presence of dipolar interactions may lead to spontaneous symmetry breaking of 
translational invariance on a scale larger than the lattice constant. 



3.2.1. Dynamical Gutzwiller approach — The time dependent version of the 
Gutzwiller wavefunction p.3p is obtained by allowing the Gutzwiller amplitudes to 
depend on time fn\t) |38| . Then the equations of motion for the amplitudes are 
readily obtained by minimizing the action of the system, given hy S — J dtC, with 

respect to the variational parameters fn\t) and their complex conjugates fn^^\t). 
The Lagrangian of the system in the quantum state j^*), is given by [40] 

C^ih ^ ' ' ' - {^\H\^), (3.6) 
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where j^*) is the time derivative of the wave function p.3p . By equating to zero the 
variation of the action with respect to /n'*\ one gets the equations 

^nin - 1) + nY,V^Af^l) ~ (3.7) 



2 



where ipi = "Pjj the sum runs over ah nearest neighbors j of site i, (hj) = 

^) is the average particle number at site j, and the total number of particles 
is given by A'^ = X^i is not difficult to verify the commutation relation 
[N, Hbb] — 0, which implies that the total number of Bosons is a conserved quantity 
for the dynamics in real time [41 \ These equations are of mean- field type, because 
for each site i the influence of the neighboring sites is taken into account in a mean 
way into the "field" <fi together with X^j^i which have to be determined self- 

consistently. Eqs. (13. 7p are a set of coupled equations, the coupling arising from the 
tunneling part, and can be written in the matrix form 

ih^f^M[f,fi,U,J]- L (3.8) 
at 



where / 



f(l) Ai) ANs) 

J jJl T ' ' ' ^ Jn 1 ' ' ' Jn„ 



is the vector of the Gutzwiller amplitudes 
ordered from site 1 to site Ns, the latter being the total number of sites. It is worth 
noticing that the matrix 7W [/, /z, U, J] is itself a functional of the coefficients / through 
the fields (pi and "Ylij^i ^ij (^i) j which have to be calculated in a self-consistent way. Let 
us clarify this point with an example. Suppose we want to solve Eq. p.8p between 
an initial time U — and a final time tf, with a given initial condition /(O). We 
discretize the time interval in N steps of size Ai, with iV finite, and define tg — sAt 
such that ts=Q = and ts=N = tf- Therefore, to calculate the solution at a certain 
point in time f{ts+i) we need to know the solution right at the preceding time J{ts), 
with which we can compute the fields that in turn determine A4[f{ts)^fJ.,U,J], and 
the solution is readily found to be 

f%^,) = e-^Minu).,.UJ]At/Hf(^^^y (3 9) 

Starting from s = 0, in + 1 steps we have determined the solution at the desired 
time tf. At the computational level, this is the simplest procedure one can implement 
to calculate the dynamics of the system. However, one needs to be careful in the 
choice of the time step At, especially for fast-oscillating dynamics. In such cases, a 
Runge-Kutta with adaptive stepsize control has proven to be more efficient. Instead in 
the case of imaginary time evolution, which requires solving stiff dynamics, the simple 
procedure described above has shown to be enough accurate and faster. 

Equations (|3.8|) can be solved in real time t or also imaginary time r = it. The 
imaginary time evolution is a standard technique that has been thoroughly used, 
because due to dissipation is supposed to converge to the ground state of the system. 
Two things are worth to be noticed. First, because the imaginary time evolution is 
not unitary, it does not conserve the norm of the Gutzwiller wavefunction, which has 
to be renormalized after each time step. Second, the total number of particles is not a 
conserved quantity anymore. For dipolar Hamiltonians the imaginary time evolution 
does not always converge to the true ground state and it gets blocked in configurations 
which are a local minimum of the energy. On the one hand this makes it a difficult task 
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to identify the ground state of such systems, and on the other hand it is a signature 
of the existence of metastable states as we wih discuss in details in Sec. H) 

3.2.2. Perturbative mean-field approach — A more convenient method to determine 
the insulating phases of a dipolar Hamiltonian is to use a mean-field approach 
perturbative in (pi. From statistical mechanics, the expectation value of the 
annihilation operator at the i-th site is given j42| by the trace 

ifii = (ai) = Tr(ajp), (3.10) 

where p = Z^^e^^^ is the density matrix operator, Z — Tr(e^^^) its normalization, 
and /3 = 1/KbT is the inverse temperature of the system. We write Hamiltonian 
(|2.54p in the form H = Hq + Hi where 

Ho ^ ^'^ni{ni - 1) - fi'^fii + '^^ninj (3.11) 

i i i=ij 

Hi=-jY,a\a„ (3.12) 

and we assume a uniform chemical potential ^. The generalization to a site-dependent 
chemical potential is straightforward. Furthermore, we assume low temperatures 
/3 — > oo, and the tunneling coefficient to be the smallest energy in the system, i.e 
J <^ U,fi^ Vij such that we can treat Hi as a small perturbation on Hq, and use the 
Dyson expansion at the first order in Hi for all the exponential operators, so that one 
obtains 

e''"'>Hie-^"'>dT . (3.13) 



We now write Hamiltonian p.l2p as a sum of single site Hamiltonians. Writing the 
annihilation operator as ai = Ai -\- ipi, we can perform the mean field decoupling on 
the tunneling term 

alttj = Alip-i + Ajipi + ipiipj + A^Aj 

~ a\Lpj + a^ipi - ipi(pj, (3.14) 

where in the last step we have assumed small fluctuations, characteristic of the Mott, 
or the deep superfluid states, and replaced AjAj ~ 0. In Hamiltonian p.l2[) we now 
replace aja_,- with the expression calculated above, we neglect terms of the order of (p'^ 
and find the mean field tunneling Hamiltonian 

i 

Given a classical distribution of atoms in a lattice such as 

i 

satisfying Hq \ $) = Eij, \ $) , let us suppose that this configuration is a local minimum 
of the energy, it can be the ground state, namely the absolute minimum, or another 
local minimum. We will be more rigorous at the end of this section regarding the 
meaning of local minimum of energy but for the moment let us refer to the common 



-P{Ho+Hi) ^ ^-I^Ho 



CONTENTS 



23 



picture of a local minimum. In the basis of the eigenfunctions of Hq, satisfying the 
relation Ho\T) — Ey\T) the partition function then takes the simple form 

Z~Tr(e-'3^«) =^(T|e-^^«|T) '-^ e"^^*, (3.17) 

IT) 

where the last limit holds because we do not trace over all the states of the basis but 
only around the state \^) , which is assumed to be a local minimum of the energy. 
Using again a Dyson expansion of the exponential of the density operator, we obtain 
the order parameter as 



~ - e'^^* / ^ Tr la,, e-t^"^)^" Hr e"^^" dr 



(3.19) 



J Tr [a, e-t^-^)^" 
rfi 

= J(p,e^^* / ^(T|a, e-f^-")^" a| e-^^^IT), (3.18) 
Jo 

which is easy to calculate. The trace is then non trivial only for |T) = |<f>) and 
|T) = -^=1$), where is integer on |$), and after the integration in the /3 i-> oo limit 
we are left with the result 

rii + 1 Tii 

E''- ^ ^ 

where the quantities E^, E'^ are defined as 

El^-^i + Un, + V^^ 

P . (3.20) 

K = M - U{n, - 1) - Vl^, 

and are respectively the energy cost for a particle (P) and hole (H) excitation on top 
of the |<I>) configuration. In the previous expressions V^^p = X^j^i ^u'^j is the dipole- 
dipole interaction that one atom placed at site i feels with the rest of the atoms in the 
lattice. We performed the integral p.lSp in the limit oi (3 i-?- oo, and in such a limit 
one finds that the integral converges only for positive values of the particle and hole 
excitation energies, namely 

U{n^ ^r) + V^<^l< Un, + V^. (3.21) 

This requirements have to be fulfilled at every site i of the lattice and they simply state 
that the configuration |<I>) is a local minimum with respect to adding and removing 
particles at any site. In the light of this statement, the restriction on the trace of 
Eq. p.l7p is now rigorous, and is in perfect agreement with the treatment done in 
[1]. Notice that if |(f>) is not a local minimum, then one finds that conditions p.2ip 
are never satisfied and the integral p.lSp indeed diverges. This treatment is of course 
also valid for |<I>), being in particular the ground state of the system. 

One finds such an equation (|3.19p . and conditions (|3.2ip for every site i of the 
lattice. The convergence conditions are simple and among them one has to choose 
the most stringent to find the boundary of the lobe at J = 0. Instead the equations 
for the order parameters are coupled due to the Cpi term, they can be written in a 
matrix form A^(/i, U, J) ■ = 0, with = {■ ■ ■ ipi ■ ■ ■) , and have a non trivial solution. 
For every /i, the smallest J for which det[A4 {fj,,U, J)] = gives the lobe of the |$) 
configuration in the J vs. /i plane. 
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3.2.3. Perturhative mean-field vs. dynamical Gutzwiller approach — The predictions 
of the perturhative mean-field treatment are in perfect agreement with the results 
of the dynamical Gutzwiller approach, since they both rely on the same mean- 
field approximation. The first looks at the stability of a given density distribution 
|<i>) = Hi of integer m atoms per site, with respect to particle and hole excitations, 
while the latter minimizes the energy of a random initial configuration with respect 
to particle and hole excitations leading to the distribution |$) if the initial condition 
is suflaciently close. However, the first method can only identify the phase boundaries 
of the insulating lobes without providing any further information on the SF phases 
outside the lobes, which can instead be explored with the imaginary time evolution. 
Nevertheless for dipolar Hamiltonians, due to the presence of many local minima of 
the energy, as we will see in the next part [13J 115] , it is very difficult to identify the 
ground state with the dynamical Gutzwiller approach. This can be achieved more 
efficiently through the perturhative mean-field approach. Therefore the two methods 
complement each other. As an example, in Fig [7] (a,b) the black lines are calculated 
with the perturhative method {Vij — 0) while the SF region outside the lobes is 
explored using imaginary time evolution showing perfect agreement between the two 
approaches. 

4. Dipolar Bosons in a 2D optical lattice 

4.1. The model 

In [331 33], we have studied the properties of dipolar Bosons in an infinite 2D 
optical lattice, mimicked by an elementary cell of finite dimensions L x L {Ns — 
sites) satisfying periodic boundary conditions. The dipoles are aligned and point 
perpendicularly out of the plane so that the dipole-dipole interaction (|2.28p between 
two particles at relative distance r becomes J7dd(r) — Cdd/(47rr'^) repulsive and 
isotropic in the 2D plane of the lattice, where Cdd is given by Eq. (|2.27|) . Furthermore 
for computational simplicity we truncate the range of the off-site interactions at a finite 
number of nearest neighbors, as shown in Fig. [HJup to the range number four. 
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Figure 9. Representation of tiie first four nearest neighbors of the site labeled 
as in the 2D lattice. 
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We have studied the phase diagram of the system described by the Hamiltonian 

77 = -J^a^aj -^//n, + ^ ^ - 1) + ^ h,hj, (4.1) 

(ij) i * e {{^3))i I ^ I 

where [/nn — Cdd/ i'^T^d^jj) is the dipole-dipole interaction between nearest neighboring 
sites, is the lattice period, and ((ij))^ represents neighbors at relative distance £ 
which is measured in units of • All other quantities were introduced previously. 



4-2. Metastability 

Let us start the discussion by introducing our definition of stability of a given classical 
distribution of atoms in the lattice, i.e. a product over single-site Fock states 

i*)=ni"^)- (4-2) 

i 

At J = 0, we define the state (|4.2p to be stable if there exists a finite interval 
= //max — fJ-min > in the fi domaiu, in which the particle (P) and hole (H) 
excitations at each site i of the lattice are positive, and the system is gapped. Using 
the dipolar Hamiltonian (|4.ip . in Sec. 13.2.21 we have calculated the particle and hole 
excitation energies of |$) to be 

where we recall V^^^ > to be the dipolar interaction experienced by one atom sitting 
at the site i of the lattice. From Eqs. (|4.3|) it is then straightforward to find Afi, if it 
exists, given by the set of inequalities 

U{n., - 1) + \^:; <^,<Un, + \^^;. (4.4) 

This is consistent with the stability conditions discussed in the seminal paper of Fisher 
et. al. [T]. Indeed, in the absence of dipolar interactions J7nn = into Eqs. (|4.4p . 
one recovers the well known conditions U {rii — 1) < n < Urii for the stability of the 
Ml(ni), with Hi particles per site. One can extend the stability analysis to small values 
of J, and for a given stable state calculate its insulating lobe with the perturbative 
mean-field approach we have developed in Sec. 13.2.21 In this context, we therefore 
define a state like (|4.2p to be metastable if it satisfies two conditions: the first is that 
the state must have an insulating lobe inside which it is gapped, and the second is 
that the energy of the state must be higher than the ground state energy. In other 
words a metastable state is a local minimum of the energy. 

In the absence of dipolar interactions C/nn = 0, no metastable states are 
found. In the low tunneling region the ground state of the system consists of Mott 
insulating lobes with integer filling factors v — Ni^/Ns (number of atoms/number 
of sites), while for large values of J the system is superfluid. In our treatment 
metastable states appear as soon as one introduces at least one nearest neighbor 
of the dipolar interaction. In fact, the imaginary time evolution, which for Bose 
Hubbard Hamiltonians with only on-site interactions converges unambiguously to the 
ground state, for the dipolar Hamiltonian (14.11) often converges to different metastable 
configurations depending on the exact initial condition. Moreover, in the real time 
evolution, their stability manifests as typical small oscillations of frequency ujq around 
the local minimum of the energy. 
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Our main results are summarized in Figs. [TU] (a,b,c), where we plot the phase 
diagram of the Hamiltonian (|4.ip for a L = 4 elementary cell satisfying periodic 
boundary conditions, for different values of the cut off range of the dipolar interactions 
respectively at one (a), two (b), and four nearest neighbors. The on-site interaction 
is given by U /J7nn = 20 and C/nn = 1 is the unit of energy. 




Figure 10. (a), (b), (c) Phase diagram with a range of the dipole-dipole 
interaction cut at the first, second, and fourth nearest neighbor, respectively. The 
thick hne is the ground state and the other lobes correspond to the metastable 
states, the same color corresponding to the same filling factor. In (c) filling 
factors range from !^ = l/8toi' = l. In the right column we present metastable 
configurations for u = 1/2 appearing at the first nearest neighbor (I), and second 
(Ila, lib), and the corresponding ground state (GS); those metastable states 
remain stable for all larger ranges of the dipole-dipole interaction. White sites 
are empty, while gray sites are occupied by one atom. 

The thick lines correspond to the ground state while the thin lines correspond to 
metastable insulating states, with the color identifying the same filling factor i/. 
The difference with the Bose Hubbard phase diagram of Fig. [71 where only integer 
filling factors v = n are present, is evident already with one nearest neighbor of the 
dipolar interaction (a). In fact the MI(1) lobe undergoes a global shift of z[/nn 
(z = 4 in the figure) towards higher values of /i, and the new fractional filling 
factor v = 1/2 appears with a ground state density distribution modulated in a 
checkerboard pattern, shown in Fig. [10] (GS). Instead, the density distribution shown 
in Fig. [TOj (I) is metastable with i/ = 1/2, and its insulating lobe is given by the 
thin line extending from 1 < /i/C/nn < 3. Remarkably, in Fig. [TUT a) the two lobes 
extending from < /i/C/nn < 1 correspond to metastable configurations at filling 
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factors 1/ — 1/4,5/16, while the two lobes between 3 < /^/C/nn < 4 correspond to 
metastable states at filling factors u — 3/4, 11/16, but no ground state is found for 
these fillings. In the region immediately outside the ground state lobes, we found 
evidences of supersolid (SS) phases, where the order parameters (pi are different than 
zero and are spatially modulated, e.g. in a CB structure. Before our work, studies of 
BH models with extended interactions have pointed out the existence of novel quantum 
phases, like the SS and checkerboard phases, but not the existence of the metastable 
states. 

Increasing further the range of the dipolar interactions leads to the appearance of 
more metastable states, as (Ila) and (lib) found at = 1/2 for two nearest neighbors 
in the range of dipolar interactions. The MI(1) undergoes a larger shift which is 
accompanied by the emergence of other insulating fractional filling factors, as shown 
in Fig. [TU] (c) , where the dipole-dipole interaction is cut at the fourth nearest neighbor 
and the ground state is a series of lobes with i/ multiple of v = 1/8. The number of 
metastable states varies depending on the parameters of the Hamiltonian and the 
filling factor; it is found to be up to 400 for C//?7nn = 20 at filling ly — 1/2, and 
up to 1500 for U/U-NN = 2 at unit filling [43]. With this picture in mind, it is now 
clear why the imaginary time evolution, which often converges to different metastable 
configurations, is very inefficient both to find the ground state of the system and to 
compute the lobe boundaries of a given metastable state. Instead, the mean-field 
perturbative approach we have derived in Sec. 13.2.21 has proven to be satisfactory for 
this purpose but it also has some limitations. In fact, all possible values of v and 
the corresponding configurations which are detectable with this method is limited by 
the size of the elementary cell. Evidently the possible filling factors of an elementary 
cell of size L are given by multiples of = 1/L^. It is worth to notice, that despite 
the inefficiency of the imaginary time evolution in finding the insulating lobes of 
the system, the corresponding equations in real time turn out to be very useful, 
for example, to compute the excitation spectrum u;(k) of the system. For a given 
metastable configuration, one can calculate w(k) from the small fluctuations 5fn\t) 

around the unperturbed metastable state coefficients /J . Writing /^'^ = /i* +^fn\t) 
into Eq. p.7p . and taking into account only linear terms in the fluctuations, one 
obtains a set of coupled equations for 6fn\t), which can be easily solved in the 
Fourier domain as explained in [44] . 

Finally, we find that there is usually a gap between the ground state and the 
lowest metastable state, which might allow to reach the ground state by ramping up 
the optical lattice under some adiabaticity condition. However, this feature is strongly 
reduced in the case of larger elementary cells because the number of metastable 
configurations and the variety of their patterns increase very rapidly with the size 
of the elementary cell L. Indeed, we have found that there exist many metastable 
configurations that differ from the ground state only by small localized defects, and 
the energy of these reduces the size of the gap. 

4-3. The lifetime 

We have studied the stability of the metastable states with a path integral formulation 
in imaginary time and a generalization of the instanton theory |45j . For any given 
initial metastable configuration I*!*) initial, we are able to estimate the time T in which 
I initial has tunnclcd completely into a different metastable state |$) final- We do this 
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in analogy with the case of a classical particle tunneling through a potential barrier 
shown in Fig. [TT] (a), with the difference that we do not have any information a 
priory on the characteristics of the potential barrier separating initial and final state. 
Nevertheless we can estimate the barrier and the time T in three steps: (i) first 
we construct the imaginary time Lagrangian of the system described by a quantum 
state I $) , (ii) we make use of a variational method on | $) with only one variational 
parameter q, and its conjugate momenta P, that interpolate continuously between 
I *&) initial and |$) final, and (iii) through the variation of q we calculate the minimal 
action 5*0, with the imaginary time Lagrangian, along the stationary path starting 
at I $) initial; this path is called an instanton path, in short instanton. It connects 
I initial and |$) final Only if the two states are degenerate, otherwise the stationary 
path connects |$) initial with an intermediate state called the bouncing point | $) bounce- 
We get an estimate of the energy barrier separating the two states by evaluating the 
Lagrangian from |$) initial to |$) final and imposing zero "momentum" P = 0, as one 
would do in the Lagrangian of a classical particle in a potential. 



Figure 11. (a) Particle in a minimum of a potential barrier, the particle oscillates 
with frequency loq around the local minimum and tunnels into the right well in 
a time T; (b) the instanton; and (c) the process for which a checkerboard state 
tunnels into the anti-checkerboard that shown complete exchange of particle with 
holes and vice versa. The process happens in a time T in analogy with (a). 



Once the minimal action 5*0 is known, the tunneling time T is readily calculated |45] 
as 



where wq is of the order of the frequency of the typical small oscillations of initial 
around the local minimum of the energy. In analogy to a classical particle tunneling 
through a barrier, the instanton has the nice interpretation of the stationary path 
connecting the two local minima in the inverted potential, as schematically represented 
in Fig. [n](b). 

In units of = 1, the imaginary time Lagrangian of a system [40j described by a 
quantum state |<I>), is given by 




final 



(4.5) 



C 



($1$) _ ($1$) 
2 



(4.6) 



with 1$) indicating the time-derivative, and H being the Hamiltonian of the system. 
In the approximation where j^) is the Gutzwiller wave function of a given metastable 
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state, we write its amplitudes as 

/W = i=(x«+zp«), (4.7) 

where Xn ^ , and pi* are real numbers which are going to be related to the variational 
parameters and their conjugate momenta in the following. For simplicity, we consider 
states with a maximum occupation number of Timax — 1 (i-c n = 0, 1), and therefore 
we have a total of iNg parameters, Ns being the total number of sites. The Lagrangian 
(|4.6p as well as the expectation value of the Hamiltonian (<i>|if|$), become functional 
of the 4:Ns parameters, namely 

1 

£[a;W,p«] Pn^^^^ + (^\^\^)- (4-8) 

z,n— 

After introducing the new coordinates d their conjugate momenta 

Pn ^ = dC/dqn^ = —ipn\ we can put the Lagrangian (|4.8p in its canonical form 

1 

-cbi^pW] = E pl'^€'^-mql:\pl% (4.9) 

where 'H[qn \ Pn'^] = — (<I'|i/|$) is a constant of the motion 0. We now want to 
reduce the dynamic described by the Lagrangian (j4.9l) to a one dimensional problem, 
described only by one variable q and its conjugate momentum P. Through the 
variation of {q, P) we want to describe the interchange between the state ji") initial and 
|^)finai, as for example the one represented in Fig. Hir e). In |Appendix A[ we show how 
to reduce the number of variational parameters to one, g, and its conjugate momentum 
P, by making use of a variational Ansatz as well as the normalization condition on 
the coefficients (|4.7p and the conservation of the total number of particles. These 
conditions are enforced through Lagrange multipliers Ac. Consequently the equations 
of motion given hj q — dV\/dP, and P — —dH/dq are governed by an Hamiltonian 
which also includes the constraints as follows 

H=n[q,P]+J2>^cCc. (4.10) 

c 

where an explicit expression for the conditions Cc is given in [Appendix A[ The action 
is then readily calculated along the stationary path of Eq. (|4.10p as follows 

So^ f C[q,P]dT= f C[q,P]^, (4.11) 
J J path 9 

with q = dH/dP from Eq. 

4-3.1. Action and tunneling time — In Fig. [12] (a, b) we plot the minimal action 
divided by the total number of sites TVs of the cell, as a function of the tunneling 
coefficient J, for two different processes. The first one (a,c), in which initial and final 
state are degenerate, shows the exchange of particles with holes in the whole lattice, 
and is sketched in the lower part of Fig. [T^] (c). There we also plot the potential 
barrier between initial and final state calculated as —H{q,P = 0). Instead, in the 



^ note that in the analogy of a classical particle in a potential V{x), the conserved quantity in the 

2m 



p2 

imaginary time would he H = ^ — — ^i^): which describes the particle's motion in the inverted 



potential. 
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second process (b,d) the final state is the ground state, i.e. deeper in energy with 
respect to the initial state, and only a few sites of the lattice exchange particles with 
holes during the process, as sketched in the upper part of Fig. [T^] (d) along with 
the potential barrier. A side remark, the point where the thick line of the barrier 
encounters the dashed line is called bouncing point | $) bounce- 




Figure 12. (a.b) Action per site and (c,d) energy barrier for the process sketched 
in panels (c,d). In both cases the initial state is configuration (lib) of Fig. 1101 
and the value J = 0.12(7nn corresponds to the tip of its insulating lobe. The first 
one (a,c) is for degenerate initial and final configurations while for the second 
one (b,d) the final configuration is energetically deeper. The difference in the two 
processes manifests also in the height of the barrier which is smaller for the second 
case, leading to a smaller action and consequently a smaller life-time. 



The action in general diverges for J — >■ indicating a divergent tunneling time T, and 
then decreases monotonically up to a minimum value in correspondence of the tip of 
the lobe, J = 0.12f/NN here, signaling a minimum life-time at the tip of the lobe, as 
expected. In between these two extreme behaviors, the action increases monotonically 
with the number of sites involved in the exchange of particles with holes; the more sites 
involved as in the case of Fig. [12] (a, c), the bigger the action is. Summarizing, from the 
figures above, we conclude that small energy differences between the initial initial 
and the final states |<i>) final and large regions of the lattice undergoing particle-hole 
exchange in the tunneling process contribute to large barriers, i.e. long life times T. 
On the contrary, for big energy differences and small regions of the lattice undergoing 
particle-hole exchange, the barrier is small. Hence, in general it is more likely for a 
given state to tunnel into a state deeper in energy, e.g., the ground state, than into 
its complementary, which implies the exchange of particles with holes in the whole 
lattice. 
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5. Multiple layers and up-down mixtures 

5.1. Dipolar Bosons in a bilayer optical lattice 

In we consider polarized dipolar particles in two decoupled 2D optical lattice layers 
(see Fig. [131), where the potential barrier between the two layers is large enough to 
prevent any inter-layer hopping. This is the simplest multi-layer structure and can be 
obtained by using anisotropic optical lattices or superlattices, which can exponentially 
suppress tunneling in one direction. 



Figure 13. Schematic representation of two 2D optical lattice layers populated 
with dipolar bosons polarized perpendicularly to the lattice plane. The particles 
feel repulsive on site U and nearest-neighbor {/nn interactions. Interlayer 
tunneling is completely suppressed, while a nearest-neighbor interlayer attractive 
interaction W is present. 



The in-plane dipolar interaction is isotropic and repulsive. The interlayer 
interaction depends on the relative position between the two dipoles, but is dominated 
by the nearest-neighbor attractive interaction W < between two atoms at the same 
lattice site in different layers. We include only nearest-neighbor (NN) in-plane (f/NN) 
and out-of-plane (W) dipolar interactions. Since tunneling is suppressed between the 
layers particles belonging to the different layers cannot mix and behave in practice like 
two different species Ltj. The problem is analogous to that of two bosonic species on a 
2D optical lattice with an inter-species attraction W < at the same lattice site, and 
intra-species repulsion C/nn- The relative strength between I/nn and W can be tuned 
by changing the spacing d± between the two layers, relative to the 2D optical lattice 
spacing d. Because of the dependence of the dipole-dipole interaction like the inverse 
cubic power of the distance, the ratio \W\/Unn can be tuned over a wide range. While 
it can be negligible for d± ^ d making the system asymptotically similar to a single 
2D lattice layer, it can also become relevant and give rise to interesting physics, not 
existing in the single layer model as pointed out in [47l l48ll49l[50l[5n[52l[53l[54] . 




d 



d 



1 



Because of this analogy we will often refer to the two layers as the two species and vice versa. 
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The system is described by the Hamihonian H = Hq + Hi , with 

E 



W^Vn«n^ (5.1) 



- J^[ala,+6l6,], (5.2) 

where a = (a, 6) indicates the two species (which in the specific case considered here 
are atoms in the lower and upper 2D optical lattice layer, respectively), U is the on-site 
energy, C/nn the intralayer nearest neighbors repulsion, W the interlayer attraction, 
J the intralayer tunneling parameter, and fi the chemical potential, as schematically 
represented in Fig. 1131 The parameters U and J are equal for the upper and lower 
layers and the chemical potentials fi are the same, since equal densities in the two 
layers are assumed. Notice that since < 0, it is necessary to have U + W > to 
avoid collapse. The symbols (ij) and indicate nearest neighbors. 

We focus on the physical situation in which the two layers are very close to one 
another {d± ^ d) such that {U + W) <C U, and small in-plane dipolar interactions 
t^NN ^ U, because in these limits particles at the same lattice site i of different layers 
pair into composites. The composites localize in a MI state for small values of the 
tunneling coefficient, while for larger values of J the pairs hop around in the optical 
lattice forming a pair-superffuid (PSF) phase [17]. Furthermore, the presence of the 
in-plane long-range interactions leads to the formation of a novel pair-supersolid phase 
(PSS), namely, a supersolid of composites [531. Finally, it is useful to introduce the 
operator of the sum of the two species number operators at each site of the lattice, 
namely 

m, = ^ \ (5.3) 

which is diagonal on a given Fock state \mi)i, with to, = {n1 + n\)/2 being the average 
occupation per layer at the site i. 



5.1.1. Low- energy subspace and effective Hamiltonian — In the limit where 
{U + W), ?7nn, J <^U, there exist a low-energy subspace spanned by the states 

\a) ^Yl\ni,n,)i, (5.4) 

i 

with equal occupation of the two species a and b at each site. These states can be 
equivalently written in terms of the pair occupations as \a) = Y[i I'nT'iji, where rui is 
the number of pairs at site i, and is given by the eigenvalue of operator (|5.3p . 

Single particle-hole excitations necessarily break one pair, and the energy cost of 
these excitations is of the order of U. Therefore, in the limit where (U + W), f/NN, J ^ 
U breaking one pair is energetically very costly. In the limit {U + W)/U — > 
0, asymptotically all states \a) become stable with respect to single-particle- hole 
excitations and develop an insulating lobe at finite J, which tend to overlap as shown in 
Fig[T5](left) by the thin blue lobes. However in this system, and for the above choice of 
parameters, single-particle-hole excitations are not the lowest lying ones. One should 
also consider two-particles-two-holes excitations, which can be accounted for within 
a low-energy theory described by an effective Hamiltonian acting on the low-energy 
subspace spanned by the states \a)-s. The validity of the effective Hamiltonian relies 
on the existence of a low-energy subspace well separated in energy from the subspace 
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of virtual excitations, to which it is coupled through the tunneling Hamiltonian Hi 
(|5.2[) via single-particle hopping. The relevant virtual subspace is obtained from the 
states I a) by breaking one composite, namely 




as schematically represented in Fig. [T31 All other states are not coupled to \a) via 
single particle hopping and hence do not contribute. 




Figure 14. Schematic representation of the low-energy subspace. The state \a) 
is uniformly occupied by one particle per site, and its energy is given by Ea- 
A second order process in the tunneling connects the two states \a), and \l3), 
through the state I7), which belongs to the virtual subspace. It is straightforward 
to notice that in the limit of {U + W), C/nn, J , the energies above satisfy the 
necessary condition Ea, Ep ^ E-^ for the existence of the subspace. 

The energy difference between the virtual states I7), and the states |a), is given by 
the sum of single-particle plus single hole-excitation energies of the states \a), which 
are of the order of C/ at J = 0, and are minimized by the width of the lobes \a) at 
finite J (see, e.g.. Fig. [15] (left)). 

Slow processes drive the system through different states of the low energy subspace 
via second order tunneling; this happens through a fast coupling with the virtual 
subspace. Since we are interested in the long time physics of the system, we have to 
average out all the fast processes and therefore we write an effective Hamiltonian in 
the subspace of pairs, and include tunneling through second order perturbation theory 
|56[ 157] . In the pair-state basis, the matrix elements of such a Hamiltonian in second 
order perturbation theory are given by 

{a\HMfi) = (a|7?o|/3) - i^(a|i?i|7)(7l^i|a) 

7 

r 1 1 1 
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where i?o, given by the interaction terms (|5.ip . is diagonal on the states |a), and the 
single-particle tunneling term Hi in Eq. (|5.2I) is treated at second order. For a given 
state I a), 

£^7,, -E„ = U+{U + W){m, - mj) + UN^Am]^^, (5.7) 

with Amjijj^ = J2{k)- — X](fc) "^fe ^ 1' where indicates the pair occupation 
number at site i as defined in Eq. (|5.3p . For (?7 + W),J7nn i7, the denominators 
E^.. — Ea are all of order U, which leads to 



rr(O) 2J^ 
er± 



i?o - — 5] [™»(m, + 1) + cjc.j , (5.8) 

(y) 

where Ci and c\ are the pair destruction and creation operators such that 
cArrii) — mAmi — 1) 

(5.9) 

cII"^^) = {rrii + l)\mi + 1). 

One can easily obtain corrections to by expanding (15. 7p at higher orders in 

(J7 + W)/U and C/nn/?^ but, as we will see, the zeroth order is already quite accurate 
to describe the physics of the system for the range of parameters we consider. 

5.1.2. Insulating lobes — We now make use of the effective Hamiltonian Hj!g derived 
above to study the ground state phase diagram of the system, starting from the 
insulating states. For every classical distribution of pairs in the lattice we can calculate 
the pair order parameters ipi — {ci) with the perturbative mean-field method derived 
in Sec. I3.2.2[ and get the expression 

2J2 \{m, + if mi 



^„ (5.10) 



where 4'i = '/'j' ^^'^ ^be energy costs of adding a pair (2P) and removing a pair 

(2H) can be calculated from the diagonal part of Eq. (|5.8p , and are respectively given 

by 

2 P 

EUJ) = + 2Um, + (2m, + 1)W + 2V^.;; - — ^(2mfe + 1) 

El^iJ) = 2^. - 2U{m, - 1) - (2m, - 1)W - 2\^:; + — Y,{2mk + 1), 

with V^^'jp — C/nn ™i being the in-plane dipolar interaction, i.e. the dipole-dipole 

interaction that one atom positioned at site i of the lattice, experiences with the rest 
of the particles belonging to the same plane. 

Using Eq. ()5.10p . one can calculate the mean- field lobes for any given 
configuration of pairs in the lattice. The ground state lobes for the checkerboard 
and doubly occupied checkerboard are shown in Fig. [15] (right) for the 0th (full lines) 
and 1st order (dashed lines) effective Hamiltonians. The comparison between the two 
shows that, for the parameters considered here, the 0th order already captures the 
physics accurately. It is worth noticing that the dependence of the energy of the 
elementary excitations is at the origin of the reentrant behavior of the lobes, which 
was predicted by exact matrix-product-state calculations for the ID geometry in [47] 
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5.1.3. Pair superfluid/supersolid — While the MI phases are predictable through 
the perturbative mean-field approach of Eqs. (|5.10p for the pair order parameters, 
to identify the SF phases, both PSF and PSS outside of the lobes, it is necessary 
to make use of the imaginary time evolution introduced in Sec. 13.2.11 based on a 
Gutzwiller Ansatz for the pair wave function. Therefore, we need to calculate the 
dynamic equations equivalent to (|3.7p in the low energy subspace described by the 
effective Hamiltonian -ff ^? of Eq. ([5 



ff 

The time dependent Gutzwiller Ansatz for the pairs, is given by 

i*)=nE/- (5-12) 

i m 

where we allow the Gutzwiller amplitudes fm (t) to depend on time. The order 
parameter is readily obtained from (I5.12p as 

= = E(™ + l)/m'^^/'mll- (5-13) 
m 

As in Sec. 13.2.11 equating to zero the variation of the action with respect to /m leads 
to the equations 



= [Um{m - 1) - 2^im + 2(Umn - — ) 



2J^\ v-^,^ , 2J 

m > {nij) zm 



J n 



2J2 



^ ^V^.m/(:^i+V^;(™ + l)/,«iJ, (5.14) 

where (m^) = I]m'™l/m^P, the fields ipi = Z](j), V'i and J^avc to be 

calculated self-consistently, and ^ = 1 the coordination number in each lattice 

layer (here z — 4). The solution of Eq. (|5.14p can be easily obtained numerically, and 
by making use of the imaginary time evolution we get the ground state phase diagram 
of Fig. [15] (right). 

To get reliable results, one should combine the Gutzwiller predictions with an 
estimate of the limits of validity of ^^ff ' beyond which the subspace of pairs looses its 
meaning. Before starting the discussion on the validity of the subspace, let us explain 
how we define the dominant classical configurations of a given state |$). It is not 
difficult to see that Eq. (15.121) can be equivalently written as 

{m} i 

where m = (mi, m^, rriNg) is a collection of the indices m at each site, and we 
have introduced the notation such that = Yl - f^} ■ The advantage of writing the 
Gutzwiller state |$) in the form (|5.15l) . lays in the product over single-site Fock states 
la) = Yli |™i)ij which is nothing but a classical distribution of atoms in the lattice. 
Therefore we can rewrite Eq. (|5.15p as 

|a>)=5].ga|«>- (5.16) 

For each point of the phase diagram, from the ground state Gutzwiller wavefunction, 
we define the dominant classical configurations with the criteria |(7^| = irii/mi"'! > 
0.02, and |/mi^P > 0.05, implying that each of the contributing f^} should also be 
sufficiently large. For each of these configurations, we calculate the lobe with respect to 
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Figure 15. (left) Pair insulating lobes for f = 0, 1/2, 1,3/2 (thick lines); Lobes 
with respect to single particle-hole excitations (thin blue lines) for the dominant 
configurations in the ground state at J = 0.05(7 and fi = —0.4375(7, namely 
nii = and rrij = 1, 2 (for i,j nearest neighboring sites). The inset shows a zoom 
of the pair phase diagram, (right) Phase diagram of the effective Hamiltonian. 
The black full lines are the semi-analytic solution of Eq. 115.101 1 indicating the 
boundaries of the insulating lobes for the checkerboard (i^ = 1/2) and the doubly 
occupied checkerboard [u = 1). The black dashed lines are the boundaries of 
the insulating lobes for 1st order expansion of HcS- The shaded area is the PSS 
phase predicted by the Gutzwiller approach. The red line indicates the estimated 
limit of validity of H^'^ (see text). The lattice parameters are C/nn = 0.025{/ and 
W = -0.95(7, which can be obtained for dj^ = 0.37d. 

single particle-hole excitationsEl- If the system at this given point of the phase diagram 
turns out to be stable against all dominant single particle-hole excitations (in other 
words, if this point is inside all selected single particle- hole lobes), H^^g is considered 
valid at that point. This procedure is shown for J = 0.Q5U and n = — 0.4375J7 in 
Fig. [15] (left), and gives the red line of Fig. [15] (right). On the right hand side of 
this red line, the low energy subspace is not well defined and therefore the effective 
Hamiltonian looses its meaning, leaving the description of the system to the domain 
of single-particle single-hole excitation theory that predicts SF and SS phases for each 
component separately. 

To summarize, in the ground state phase diagram of Fig. [15] (right) we identify 
four different types of phases characterized by different values of the order parameters, 
the single-particle order parameters for each species (pf — (a;), ip^ = {bi), and the pair 
order parameter -0i — (ci). 

(i) The Mott insulating checkerboard phases (MI) characterized by fractional filling 
factors I' and vanishing order parameters cpf = ^p^ ^ ^pi — inside the lobes. 

(ii) The pair-superfluid phase (PSF) in which the single particle order parameters are 
zero (ff = ip^ = 0, while the pair order parameter is uniformly non-zero ipi 0, 
signaling a finite fraction of superfluid density of the pairs. 

(iii) The pair-supersolid (PSS) characterized by vanishing single-particle order 
parameters (pf = (p^ — 0, and non vanishing pair order parameter ipi ^ 0, 
coexisting with broken translational symmetry, namely, a modulation of both 
density and order parameter on a scale larger than the one of the lattice spacing, 
analogously to the supersolid phase. 

* We have checked that the validity region is not strongly modified upon small changes in these 
conditions. 
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(iv) We infer superfluid and supersolid phases of both components SF^ and SF^ (SS^ 
and SSb respectively) with non vanishing single species order parameters ip'^ ^ 
and (^^ 7^ 0, as well as non vanishing ijji 0, on the right hand side of the red line 
of Fig. [15] (right), which we estimate to be the limit of validity of the low energy 
subspace for the parameters considered here. 



5.2. Up-down mixture in a 2D optical lattice 

In most cases, if not all cases considered so far, dipolar gases are polarized, i.e. 
magnetic or electric dipoles point in the same direction. This, however, does not 
have to be always the case. A prominent example concerns molecules that follow 
the Hund's rule (a) [58]. For such molecules the electric dipole is either parallel or 
anti-parallel to the direction of the magnetic moment. Thus, if one takes a sample of 
such molecules and polarizes it magnetically, say in the up direction, one will obtain 
in this way a dipolar gas with certain fraction of electric dipoles pointing up, and the 
remaining fraction pointing down. In fact there was recently a spectacular progress 
in cooling and trapping of magnetically confined neutral OH molecules. This progress 
opens important perspectives to study novel quantum phases with ultracold dipoles 

[Mlleo]. 

In [61] we consider a sample of dipoles in the presence of a 2D optical lattice, 
and an extra confinement in the perpendicular direction. The dipoles are free to point 
in both directions normal to the lattice plane, which results in a nearest neighbor 
interaction either repulsive for aligned dipoles, or attractive for anti-aligned one, as 
shown in Fig. 1161 The system is described by the Hamiltonian 



H = 



E 



-TT^T-t [ri., - 1) + -rrrii n. 



Hern, 



E 



Unn 



-^E 



(5.17) 

where a = (a, 6) indicates the two species, i.e. dipoles pointing in the up and down 
direction perpendicular to the 2D plane of the lattice, respectively. Uaa and Ubb are 
the on-site energies for particles of the same species, while Uab is the on-site energy for 
different species. The long-range dipolar interaction potential decays as the inverse 
cubic power of the relative distance r^ , which we express it in units of the lattice 
spacing. For computational reasons, in most theoretical approaches the range is cutoff 
at a certain neighbors. In the present work, we consider a range of interaction up to 
the 4th nearest neighbor. The first nearest neighbor dipolar interaction is attractive 
(repulsive) for particles of the same (different) species, with strength Unn > (or 
respectively —Unn), and we consider an equal tunneling coefficient (J) for both 
species, while the densities of the dipoles pointing upwards and downwards are fixed 
by the corresponding chemical potentials /io-. 



The on-site interactions have two contributions (see Eq. I2.55P : one is arising from 
the s-wave scattering length, and the second one is due to the on-site dipole-dipole 
interaction. We consider that the s-wave scattering length is independent on the 
orientation of the dipoles. Instead, the on-site dipolar contribution U^d depends both 
on the orientation of the dipoles and on the geometry of the trapping potential, and it 
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Figure 16. Schematic representation of a 2D optical lattice populated with 
dipolar Bosons polarized in both directions perpendicular to the lattice plane. 
The particles feel repulsive intra-species Uaa, Utb, and inter-species Uab repulsive 
on-site energies. The nearest-neighbor interaction is repulsive C/nn > for aligned 
dipoles, while it is attractive — C/nn for anti-aligned particles. The hopping term 
J is equal for both species. 



can be varied by changing the ratio between the vertical to the axial confinement. For 
simplicity we will focus on the specific case of a spherically symmetric confinement, 
where the on-site dipolar interactions average out to zero Udd — 0, and the resulting on- 
site interactions are all equal to U. We consider the case of dipole-dipole interactions 
to be 600 times weaker with respect to the on-site interaction, i.e. C/nn = [7/6000. 

The properties of the system are conveniently extracted using the operators given 
by the sum (filling factor) and by the difference (imbalance) of the two species number 
operators at each site of the lattice, namely by 

nf + fi^i nf — , , 

= m, = (5.18) 

which are simultaneously diagonal on a given Fock state m) ;. Notice that the 
eigenvalues of these two operators are not independent. In fact, by fixing v the 
eigenvalues of rhi can only assume 2v ~\-\ values given by m = {~^^ 1, in 
complete analogy with the angular momentum operator Sf and its projection along 
the z axis S^ . It is also useful to introduce the average magnetization of the system, 
defined as 

i 

where Ng is the total number of lattice sites. 

Substituting Eqs. (jSlS]) into Eq. ([STZ)) 
Hamiltonian as H = + HJ^ + H'("', where 



allows us to express the system 



K = 



2Uv,, 



(5.20) 



tl In standard experiments with ^^Cr, which features a magnetic moment of ^ = this ratio 

is given by (7nn — f7/400, for an optical lattice depth of 20-Er, where i?R is the recoil energy at 
A = 500nm. 
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Tjm 
^0 



E[-2M--. + 2C/^.^.5:!^] (5.21) 

Hr=-jJ2[a\a,+brb,]. (5.22) 

iv) 

In Eqs. (|5.20l5.2ip . we have introduced the chemical potentials 

M± = ^ , (5.23) 

which respectively fix the eigenvalues of the filling factor (+), and the imbalance 
operators (— ) in Eq. (|5.18|) . In the following we consider Hi™ to be a small 
perturbation on the interaction terms (15.201) and (|5.2ip . 

5.2.1. Low-energy subspace and effective Hamiltonian — In the limit where U ^ 
Unn-, there exist a low-energy subspace spanned by the states 

\a) =l[W,m,),, (5.24) 

i 

with uniform total on-site occupation 2i/. Single particle hopping changes the total 
on-site population, the energy cost of these excitations is of the order of the on-site 
interaction energy U, and becomes very large in the limit where U (Unn, J)- Thus, 
a successful description of such a system is obtained through an effective Hamiltonian 
Heff restricted to the low-energy subspace, where single-particle hopping is suppressed 
and tunneling is included at second order in perturbation theory. 

This situation is in fact similar to the one discussed in Sec. l5.1l for a bilayer optical 
lattice, and therefore we apply the same technique to compute H^s- In the basis of 
constant total on-site population 21/, using Eq. (|5.6p we calculate the matrix elements 
of such a Hamiltonian at second order in perturbation theory, where the subspace of 
virtual excitations I7) is obtained from the states |a) via single particle hopping, as 
schematically represented in Fig. [iTl 

lv> 



VWV 



WW WW 



Figure 17. Schematic representation of a two-particle hopping between the states 
q) and 1/3), belonging to the low-energy subspace at !^ = 1/2. These states are 
coupled through virtual excitations to the states I7) by single-particle hopping. 

For a given state \a), 

E^^^ ~E^^U + UnnAtu^j,, (5.25) 

with AtoJ;Jj^ = I]fe^j 2r7ifc/|riA:p - J2k ^i 2m fc/|rjfcp - 1, where indicates the 
population imbalance at site i of Eq. (|5.18p . For U ^ Unn, the denominators 
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Ej.^ — Ea are all of order [/, which leads to 



2 j2 p 1 

^o"-— Ei™J- + ^Ic,J, (5.26) 

where Ci = a^SJ and cj = d|6i are composite operators, corresponding to the creation 
of a particle of one species and a hole of the other species, such that 



Ci\iy,mi) = + 1) - miinii - l)\iy,mi - 1) 

cj|z/, mi) = \/v(i> + 1) - mj(TOi + \)\v, rUi + 1). 

5.2.2. Insulating lobes — In the limit where U ^ ?7nn, asymptotically all classical 
states I a) become stable with respect to single-particle-hole excitations and develop 
an insulating lobe at finite J. As single particle hopping changes the total on-site 
population, it breaks the translational invariance of the ground state with respect to 
the total on-site occupation 2z/. The energy cost of these excitations is of the order 
of the on-site interaction energy U, and is therefore very large in the limit where 
?7» {Unn,J). 

Instead, the lowest-lying excitations within the subspace are obtained by adding 
(PH) or removing (HP) one composite, made of a particle of the upper-polarized 
dipoles (species a) and a hole of the lower-polarized dipoles (species b), at the i-th site 
of the lattice. For any given configuration |a), one can calculate the corresponding 
energy costs from the diagonal terms of the effective Hamiltonian (j5.26p . which are 
respectively given by 



2 



k^t {k}i 



i;^p(J) = 2/i_-4C/NNE 



JTlfc 4J^ 



(5.28) 



fc/i {k)i 



^ mk. 



The order parameters ipi = (c^) for \a), are readily calculated as in Sec. I3.2.2[ and 
satisfy the equations 



, 2.P 



h'{v + 1) — mi{mi + 1) v{v + 1) — mi{mi — 1) 



where ipi = ^j- With Eqs. (|5.29l) one can calculate the mean- field lobes of any 

distribution of atoms in the lattice la). 

In Fig. [T8]we plot the ground state insulating lobes calculated in this way for i/ = 1/2 
(left) and i/ = 1 (right). For all filling factors i/, we find an anti- ferromagnetic (AM) 
ground state (i/, M = 0). The AM insulating lobes are symmetric with respect to the 
— axis, and present a spatial distribution of alternating sites occupied by particles 
of species a and b resembling a checkerboard structure. Remarkably, the larger the v, 
the more stable is the AM ordering with respect to hipping the direction of a dipole. 
By increasing the absolute value of fi_ we find RM ground states with rational values 
of the average magnetization, corresponding to M = (±2z/, ±4j/, ±6z^)/8. The exact 
fractional values of M in the ground state, depend on the cutoff range in the dipolar 
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Figure 18. Ground state of a 4 X 4 square lattice satisfying periodic boundary 
conditions, for u = 1/2 at fi^ = 0.5U (left), and f = 1 at /i^ = lAU (right), and 
f^NN = U/600. The text in parentheses (v,M) indicate the filling factor v, and 
the average magnetization M, respectively. 



interactions, and on the size of the lattice. We have used a 4 x 4 elementary cell with 
periodic boundary conditions, and dipolar interaction range cut at the 4th nearest 
neighbor. By considering more neighbors in the interactions, and larger lattices, we 
expect to find RM states appearing at all rational M, asymptotically approaching a 
Devil's staircase as recently shown in [79l [13]. Finally we find a FM ground state 
(z/, M = ±z/), in which only particles of one type are present. 

It is worth noticing that the insulating lobes calculated in this way, do not contain 
any dependence on which does not enter into Eqs. (j5.29p . Therefore, for any given 
value of /u+j in order to obtain the ground state phase diagram one has to compare 
the energies of the ground state configurations at different v. Using the effective 
Hamiltonian (|5.26|) . for any value of J, and /x„, we compare the energies of the 
ground state configurations for different and select the state with the smaller energy. 
In this way we have obtained the phase diagram at J = shown in Fig. 1191 



5.2.3. Counterflow superfluid/ super solid — In the low-energy subspace at constant 
V, the Gutzwiller Ansatz on the wave function of the system reads 

V 

l'i')=n E /^^l'^'™)- (5-30) 

i m— — 

where we allow the Gutzwiller amplitudes f^%(t) to depend on time. We obtain the 
equations of motion for the amplitudes by minimizing the action of the system with 
respect to the variational parameters f^!m{t), 



0). 
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Figure 19. Ground state of the system at J = 0, calculated for a 4x4 elementary 
cell satisfying periodic boundary conditions, and (7nn = U/600. The text in 
parentheses {u, M) indicates the filling factor u and the average magnetization 
M. 



U 



i,* VK'^ + l)-m(m + l) fl%^,\ , (5.31) 



where (m^) = I]m=-i/ ^^'^ ^elds = V'j' 



P have to be calculated in a self consistent way, and the order parameter 
tpi = (<I>|ci|<I>) is given by 

V'.= E v/K'^ + l)-™(m + !)/;« (5.32) 

m— — 1> 

We solve Eqs. (I5.3ip in imaginary time t = it, and in Fig. [TS] we show the ground 
state phase diagram of the system for = 1/2 (left) computed in this way. 

For V = 1/2, in the region immediately outside the insulating AM lobe, depending 
on the values of J and we find either super-counter-fluid SCF or counterfiow- 
supersolid CSS. In the SCF phase, while the single-particle order parameters vanish 
(cii) = {hi) = 0, the composite order parameters are non-zero (ci) 7^ 0, indicating the 
presence of counterflow [57j. The CSS is characterized by vanishing single-particle 
order parameters {hi) = {hi) = 0, and non-vanishing composite order parameters 
{ci) 7^ 0, coexisting with broken translational symmetry, namely, a modulation of 
both mi, and (c^) on a scale larger than the one of the lattice spacing, analogously to 
the supersolid phase. Note that we don't find any evidence of CSS at ^_ — 0, which 
indicates that a finite imbalance between the two species is a necessary condition for 
the system in order to sustain CSS. Finally, with a similar method described in Sec. 
15.1.31 we estimate the limits of validity of the effective Hamiltonian to be given by the 
vertical thick lines of Fig. [TBI 
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6. Path Integral Monte Carlo and the Worm algorithm 

Any physical system consisting of N non-relativistic particles can be in principle 
described by the many-body Schrodinger equation. In three dimensions (3D), the 
number of degrees of freedom in the Schrodinger equation becomes 3 times N. For 
typical physical systems such as electrons in conducting materials or BEC that have 
a large number of constituents N, the Schrodinger equation becomes difficult to solve 
exactly in a reasonable amount of time even for parallel computing. Monte Carlo 
methods overcome this problem, they allow for a description of the many-body system 
relying on repeated random sampling, at the cost of statistical uncertainty which can 
be reduced with more simulation time. The typical basic steps of a Monte Carlo 
algorithm can be summarized as follows 

(i) Define the configuration space (here by configuration we mean a collection of 
indices, see below). 

(ii) Generate configurations randomly and accept them with a certain probability 
which depends on the specific problem. This is called the updating procedure. 

(iii) Perform a computation, i.e. calculate quantities of interest, based on the 
randomly generated configurations. 

(iv) Add the result of the computation to the final result. 

In the spirit of the Metropolis-Hastings algorithm [67l |68] , two requirements must 
be satisfied: a) ergodicity, that is, given an initial configuration it has to be possible to 
reach any other allowed configuration via the updating procedure; b) the probability of 
having a certain configuration appearing in sums which calculate quantities of interest 
has to be proportional to its Boltzmann weight. 

There is a large class of Quantum Monte Carlo methods that can simulate 
quantum many-body systems, like for example the Variational Monte Carlo ^62', '63], 
the Diffusion Monte Carlo [H EH [65] , the Path Integral Monte Carlo [66, .69, 70], 
auxiliary field Monte Carlo [7TJ [72] , etc. Most methods aim at computing the ground- 
state wavefunction of the system, with the exception of Path Integral Monte Carlo, 
and finite-temperature auxiliary field Monte Carlo, which calculate the density matrix. 
The results presented in this work are based on the Path Integral Monte Carlo 
(PIMC) and the Worm algorithm (WA), which was originally developed by Prokof'ev, 
Svistunov and Tupitsyn (55] [7D] . 



6.1. Path Integral Monte Carlo 



Consider a system described by the Hamiltonian H = Hq + Hi , where Hq is diagonal 
in the basis of eigenstates {la}} satisfying the eigenvalue equation 

Ho\a)=E^\a), (6.1) 

and Hi is non-diagonal. The thermodynamic properties of the system at equilibrium, 
can be derived from the partition function which is given by the trace of the density 



matrix operator Z = Ti 



, where /3 = 1/KbT is the inverse temperature and 



Kg the Boltzmann constant. In the interaction picture |42| one may write 



Z = Tr 



,-/3(Ho+Hi) 



= Tr 



(6.2) 



where 7^ is the time-ordering operator, Hi(t) = e'^ "Hie " is the non-diagonal 
part of the Hamiltonian expressed in the interaction picture, and the variable r is 
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usually called the imaginary time^^. One can write the partition function using the 
Feynman path integral formulation, and by Taylor expanding the second exponent in 
the right-hand-side of Eq. (j6.2p one gets 

Z= Ve-'3^°(a|i- / dTHi{T) + 

+ dT„, .. / dTiHi{T^) ..Hi{Ti)\a), (6.3) 

m=2 

where the integrals are ordered in time and the sum over the states |a) comes from 
the trace. Now we explicitly make use of the completeness property of the {|q;)} basis, 
and insert m — 1 identity operators 1 = J2a between the products of ffi(rm) 

operators, therefore we can write 

(a|77i(T™)..i?i(ri)|a)= Yl Hr-\r^) ■■ H'^'"^ {T2)H'^'''{n), (6.4) 

ai ,..,am— 1 

where the matrix elements 

if«'"(r) = e^^-'Hf^e-^^" = {a'\Hi\a)e-^^^''-^-'\ (6.5) 

contain both diagonal (Ea) and off-diagonal (i?f ") matrix elements. We now insert 
the last equation into expression (|6.3p and get the final expression for the partition 
function 

Z= Ve-'^^Wl- / dTiJrM+ (6.6) 

+ ^(-1)" / dr.„.. / dn i?r"-^(r.„)-j^r"(n)}, 

, 9 "'0 Jo „ , 



ai ,. . ,Qm— 1 



which contains only matrix elements of the operators Hq and Hi. Therefore, by using 
this formalism of path integrals, the calculation of the partition function reduces to a 
classical problem since only scalars enter into Eq. (|6.6I) . but we have payed the price 
of the extra dimension r. In other words, the original d-dimcnsional quantum system 
is equivalent to a (d -f l)-diniensional classical system. 

It is worth noticing that since the partition function is a trace, periodic boundary 
conditions in the imaginary time r must apply. This is easily understood by looking 
at the m-th order term of Z, which contains the product of m matrix elements 
iJ"""~^(T,„) .. _ff"^"(ri) that are ordered in time from the first at ri, to the last 
at Tm- Therefore, for any given a in the trace, the first matrix element brings a to 
some ai at time ri > 0, while the last matrix element brings am-i back to a at time 
TVn < /3- All the possible configurations which are periodic in imaginary time and that 
enter into the expression for the partition function Eq. (16. 6p . define the configuration 
space spanned by a PIMC algorithm. 



6.1.1. Path Integral Monte Carlo and the 2D extended Bose-Hubbard model — We 
now consider a 2D system oi LxL sites filled with polarized dipolar Bosons, we assume 
spatial periodic boundary conditions and the dipoles to be polarized perpendicularly 
to the 2D plane as explained in Sec. The system is therefore described by the 

ft This is because by replacing t = it, with t being the real time, the operator e~^^ becomes the 
usual time-evolution operator in quantum mechanics. 
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extended Bose-Hubbard Hamiltonian (|4.ip . which, to be consistent with the notations 
in our pubhcation [79] . we rewrite in this form 



H=-jj2brb,+j2 



u 

—ni(ni - 1) - ^iJii 



where bl (bi) is the boson creation (annihilation) operator at site i, fii = b\bi is the 
number operator, V = D/a^ > Q \s the dipole-dipole interaction strength D divided 
by the lattice spacing a, = |j — j| is the distance between two sites of the lattice, 
and = ^ — fli^ contains the chemical potential /i which fixes the number of particles, 
and the curvature Q of an external harmonic confinement. 

We choose to work in the basis of the interaction term of the Hamiltonian (|6.7p . 

i.e. Fock states \a) = Yli \n'i)i of localized particles in the Lx L square lattice, where 
Hi is the occupation number at site i. Therefore in this basis, the diagonal matrix 
elements entering Eq. (j6.5p take the form 

= I ^ ruin, f^^^^ + ^ E 7?^- (6-8) 

i i i<j '^J 

The off-diagonal ones are given by the expression 

- fff'" - J(a'|6l6,|a) = J^{nf + l)nf, (6.9) 
and they connect states jo;') and \a) that differ only in the occupation number of the 

b - 

two nearest neighboring sites i and j, namely \a') = —j=^===\a) with nf being the 

integer number of particles at the i-th. site of the state \a). 

To write the partition function for the 2D extended Bose-Hubbard model, we 
notice that the first order term vanishes since the matrix elements (|6.9p are off- 
diagonal, i.e. H"" — 0, and due to the geometry of the system (2D square lattice) 
it is not difficult to see that all the terms with an odd value of m also vanish, owing 
to periodic boundary conditions in imaginary time. Therefore by rearranging the 
exponentials and renaming a = ao, we get to the expression 

oo 

^.BH = E + E ^"^^ ^ (6.10) 

ao m— 2 

X / dr^ .. / dri ^ exp I- l3Ea„ - ^ Ea^Tp+i - Tp)j, 

where is a product of m square root factors coming from Eq. (j6.9p and we have 
introduced tq = t„i to compact the notation. We can compact further the notation 
by defining Am=o — 1, and keeping in mind that for m = there is no summation in 
the exponent of Eq. (|6.10p we then write the partition function as follows 

oo 

= E E ^ (6.11) 

m— ao,ai,..,(im-i 

/ drm .. / dri exp \ - PEa^ - Ea^{Tp+i -Tj^\- 
Jo Jo ^ J 

From the last expression one can formally write 

^cBH^E^-' (6.12) 
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with Wi, being the weight of each configuration v = [m,ao{Ti),ai{T2), ■.,atn~i{T„i)]. 
A configuration can be pictoriahy represented as in Fig. 1201 with imaginary time t 
on the horizontal axis, and lattice sites on the vertical axis. 




I ■ >- 



Figure 20. Schematic representation of one configuration which enters the 
calculation of Zobh- Each line is called a worldline and it represents a trajectory 
of a particle in imaginary time. Configurations have to fulfill periodic boundary 
conditions in the imaginary time r, owing to the definition of the partition function 
as a trace. Vertical arrows correspond to changes of the state of the system, and 
are called kinks. In the sketch, the thickness of intervals between the kinks shows 
the number of particles: the dashed black line is for n particles while the solid and 
bold blue lines have occupation numbers equal to n + 1 and n + 2 respectively. 



Each line represents a trajectory of a particle in imaginary time and is called a 
worldline. The latter has to close on itself owing to the fact that the partition 
function is a trace. Moreover if one assumes spatial periodic boundary conditions one 
can imagine the configuration of Fig. [20] to be wrapped on a torus (in the case of one 
dimensional systems) . We call the phase space of all possible configurations the closed 
path configuration space (CP). 

If we cut one configuration at a certain instant in imaginary time, we get the 
system in a particular quantum state. The points in imaginary time where the system 
changes state are called kinks, which in Fig. [20] are represented by vertical arrows. 
A configuration with a number of kinks equal to m, contributes to the m-th order 
term of the partition function Eq. (|6.1ip . and it is straightforward to see that there 
exist an infinite number of different configurations with the same number of kinks, the 
difference being the time at which the kinks take place and/or the different states they 
connect. The updating procedure of a PIMC algorithm therefore consists of changing 
the number of kinks and/or their position in imaginary time. We will discuss the 
updating procedure specifically for the Worm algorithm in the next section. 
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6.2. The Worm algorithm 

The Worm algorithm, originally developed by Prokof'ev, Svistunov and Tupitsyn 
[Bni[7D], works in an enlarged configuration space, in which one allows one disconnected 
worldline, the worm, drawn as a red line in Fig. 1211 This is equivalent to work in 
the Grand-Canonical ensemble, as we shall discuss in Sec. 16.2.11 The disconnected 
worldline allows to efficiently collect statistics for calculating the Matsubara Green 
function, defined as 

G(j,r) = (r.Si+j(To + r)6j(ro)), (6.13) 

where T5- is the time-ordering operator, tq and r are two points in imaginary time, i 
and j are two sites of the lattice, and the symbol (.) stands for the statistical average of 
the expectation value of an operator. Due to space and imaginary time translational 
invariance of the system, the Green function Eq. (|6.13p does not depend on i and 
Tq. The configuration space of the Matsubara Green function is called the CPg space, 
and it is easy to see that the only difference between configurations contributing to 
the partition function .^obh and those contributing to the Green function G is that, 
for the latter, one of the worldlines starts at (i, tq) and ends at (j + j, tq + r), i.e. the 
worldline is disconnected. 




Figure 21. Configuration of the CPg space, tlie red disconnected line represents 
the worm. 



6.2.1. Updating procedures — Let us now discuss the updating procedure of the 
Worm algorithm, that is when the system is in a certain configuration v and the 
algorithm has to generate randomly a new configuration v' to collect statistics for 
evaluating the observables of interest. 

Apart from the creation of a worm, which is done in the CP space, all other 
updates are done in the CPg space through the two ends of the worm. Nearly all 
updates are done in pairs. One can picture the updating scheme as sequence of 
'drawing' and 'erasing' procedures, happening at the end points of the worm. Below 
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we list and describe the four types of updates the Worm algorithm uses to simulate 
single component Bose-Hubbard models. 

Creation of a worm — Creating a worm is the only update performed in the 
CP space, therefore the starting point is a configuration i' belonging to CP. Each 
configuration can be thought as divided into intervals, each interval being delimited 
by kinks shown in Fig. [22] as crosses (or worm extremities for configurations in the 
CPg space). In the create worm update one of the intervals of v, delimited by r^-m and 
Tinax (see interval ni in Fig. B^ . is randomly selected. Then the algorithm suggests at 
random two points ti and T2 within ni, which will be the worm extremities (indicated 
by plain dots in Fig. I22p . with the constraint Tmin < ti < T2 < t^hx- With equal 
probability one suggests to draw or delete the piece of worldline delimited by ti and 
T2, with the constraints that the resulting configuration belongs to the Hilbert space, 
i.e. it is not possible to erase from an empty interval or to draw on an interval which 
has reached the maximum occupation number allowed, if any. The worm is therefore 
created and all other updates will take place through its two extremities. 




n T2 n T2 

)C ■ ■ )( 

6 St St s 

Figure 22. Creation/Deletion of a worm. In the create worm update an interval 
is randomly selected (top), and two points ti and T2, which will become the two 
extremities of the worm, are randomly chosen. Then one can either delete a piece 
of worldline (bottom left) or draw a piece of worldline (bottom right) with the 
constraint that the Hilbert space is observed. 



Deletion of a worm — In analogy, the update opposite to creation of a worm is 
the deletion of a worm. It can only take place in the CPg space and only if the two 
extremities of the worm belong to the same interval. 

Time shift — This is the simplest of the updates and it consists of moving one of 
the extremities of the worm in the imaginary time direction, such as to lengthen or 
shorten the size of the worm. The algorithm selects at random the imaginary time 
instant to which the extremity of the worm will be moved. 

Space shift — This update changes the number of kinks and it consists of creating 
or deleting a kink to the left (space shift left) or to the right (space shift right) of 
the worm extremity. Fig. I23f a) shows the creation/deletion of a kink backward in 
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imaginary time, i.e. the space shift left. In the creation update a nearest neighbor of 
the site to which the worm extremity belongs is selected at random and the kink is 
inserted at an imaginary time instant within the interval delimited by T^^in E^nd r,nax 
[see Fig. I^HTa)]. i.e. with the requirement that the created (or deleted) kink does not 
interfere with any other interval. 



[a) 
< ► 

■'"max 
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Figure 23. Sketch of space shift updates that create or delete kinks. In the 
space shift left (a), one kink is created/deleted backward in the imaginary time, 
i.e. to the left of the worm extremity; in the space shift right (b) the kink is 
created/deleted to the right, i.e. forward in the imaginary time. 



The last update, the space shift right shown in Fig. [SBT b). is equivalent to the left one 
with the only difference that the kink is inserted or deleted to the right of the worm 
extremity, i.e. forward in imaginary time. 

These are all the updates performed by the WA. It is straightforward to see 
that the WA works in the Grand Canonical ensemble, i.e. allows to change the 
number of worldlines present in the configurations. The chemical potential becomes 
an input parameter which fixes the average particle number. For example, suppose the 
algorithm starts with an initial configuration v of zero particles in the system. From 
this configuration the only possible update is to create a worm by drawing a piece 
of worldline. Then, trough the space shift and time shift updates, the worldline will 
eventually close on itself, corresponding to the insertion of one particle in the system. 

6.2.2. Advantages of the Worm algorithm — The updates described above are 
all local and allow to draw/erase any line, and create kinks between the sites. 
Although only configurations belonging to the CP space contribute to the evaluation 
of the partition function, by using the enlarged configuration space CP + CPg the 
intermediate configurations with one disconnected loop allow to efficiently collect 
statistics for the Green function. For an algorithm working in the CP space only, 
instead, collecting statistics for the Green function results computationally very 
expensive. 
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Another advantage of the WA is that it does not suffer from critical slowing 
down in the vicinity of a critical point. In the critical region, a system develops long 
range correlations, and in most cases an algorithm based on local updates results very 
inefficient in simulating such a system for which the relevant degrees of freedom are 
non-local. This results in the divergence of the autocorrelation time with the system 
size. Although the WA performs local updates, it overcomes this problem by using the 
drawing and erasing updating procedures through the worm ends, which are directly 
linked to the critical modes (long range order in G{j,T)). As a result generating 
independent configurations in the critical region is very efficient. 

The WA is also efficient in sampling topologically different configurations and 
configurations which are separated by an energy barrier. This property is a necessary 
condition in order to maintain ergodicity. An example of two topologically different 
configurations is shown in Fig. I24[ where a one-dimensional system with one particle 
(worldline) is considered. Periodic boundary conditions in time and space apply, i.e. 
the system is a torus where the bottom and top facets of the cylinder are glued 
together. Fig. HDJa) represents a configuration with zero winding numbers, i.e. the 
worldline does not 'wind' in space. Fig. I24f d). instead, represents a configuration 
with one winding number, i.e. the worldline winds once in space. An algorithm 
based on local updates which only works in the CP space would not allow to sample 
configurations with different winding numbers, unless a global update which introduces 
a winding number at once, is introduced. The WA, instead, can easily go from 
configuration of Fig.lMTa) to configuration Fig. \24[d) (see a sketch in Fig. [24lb)-(c')'). 

Being able to sample configurations with different winding numbers is crucial in 
order to simulate SF systems. It was shown in [T^, that the superfluid stiffness can 
be extracted from the statistics of winding numbers 

Ps 



TiW 
'dl? 



(6.14) 



where T is the temperature, L the system size, d the dimensionality, and = 
SiLi ) with Wi being the winding number in the coordinate i. 
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Figure 24. Onc-dimcnsional system with (a) zero and (d) one winding 
number(s). {b)-(c) sketch of how the WA is able to go from (a) to (d). 
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Figure 25. Phase diagram corresponding to the Hamiltonian Eq. 1)6. 7|l as a 
function of /i and J at zero temperature. Lobes: Mott solids (densities indicated); 
SS: supersolid phase; SF: superfiuid phase. DS: parameter region where devil' s 
staircase is observed. Panels (b-d): sketches of the groundstate configuration for 
the Mott solids in panel (a), with p = 1/2, 1/3 and 1/4, respectively. 



7. Quantum Monte Carlo studies of dipolar gases 

Quantum Monte Carlo is one of the most powerful methods we have to study 
equihbrium properties of strongly interacting many-body quantum systems. In the 
literature, there is a large amount of work devoted to the study of dipolar gases with 
Quantum Monte Carlo techniques. From self-assembled floating lattices, provided by 
trapped polar molecules [74j, to the possibility of tuning, and shaping the long-range 
interaction potential of polar molecules [7S], to self-organized mesoscopic structures 
of matter waves in zigzag chains [65] . to the spectrum of the elementary excitation 
that can exhibit a roton minimum |76i 177] . to the emergence of an emulsion phase 
in triangular lattices [78]. The ones listed above are just a few of the outstanding 
properties of dipolar gases, which have been investigated with various Monte Carlo 
techniques. 

Based on the Path Integral Monte Carlo and the Worm algorithm, in [TS], we 
have studied the ground state properties of dipolar hard-core Bosons confined in a 2D 
square lattice of linear size L, satisfying periodic boundary conditions. The system is 
described by the extended Bose- Hubbard Hamiltonian (|6.7I) . where no cut-off in the 
dipolar interaction potential is used. 
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7.1. Incompressible and supersolid phases 

The incompressible and supersolid phases are both characterized by a finite value of 
the structure factor, defined as 

r,r' 

with k the reciprocal lattice vector, the density at position r, and N the total 
number of particles. While for the incompressible phases the superfiuid fraction 
vanishes Ps ~ 0, the supersolid phase is characterized by a finite value of Ps, indicating 
the presence of superfiuid. 

Our main results in the absence of harmonic confinement = 0, are summarized 
in Fig. [25l where we show the zero temperature phase diagram of the system, in the 
J vs. fj. plane, in the range J/V > 0.02, and 1 < fi/V < 6 indicated by the unshaded 
area. For finite J, three main solid Mott lobes emerge with filling factor p — 1/2, 1/3, 
and 1/4, named checkerboard (CB), stripe (ST), and star (SR) solids, respectively. 
The corresponding groundstate configurations are sketched in panels (b-d). We find 
that the CB solid is the most robust against hopping and doping, and thus it extends 
furthest in the J vs. /i plane. 

For large enough J/V, the low-energy phase is superfiuid (SF), for all p. At 
intermediate values of J/V, however, we find that by doping the Mott solids either 
with vacancies (removing particles) or interstitials (adding extra particles) a supersolid 
phase (SS) can be stabilized, with coexisting superfiuid and crystalline orders. Instead, 
we find no evidence of SS in the absence of doping. The green shaded area above and 
below the CB lobe boundaries in Fig. [25l correspond to a SS obtained by doping the 
CB crystal with interstitials, and vacancies respectively. Remarkably, the long-range 
interactions stabilizes the supersolid in a wide range of parameters, in fact for one, or 
two nearest neighbors in the dipolar interaction range, no stable CB SS was found for 
P < 1/2 [H. 

Interestingly, we find evidence for incompressible phases in addition to those 
corresponding to the lobes in Fig. [25l This is shown in Fig. [26l where the particle 
density p is plotted as a function of the chemical potential p. 
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Figure 26. p vs. fi. (a): Solids and SS for a system with linear size L = 12 
and J/V = 0.05. Some p are indicated, (b): SF and vacancy-SS for L = 16 and 
J/V = 0.1. 



CONTENTS 



53 



In Fig. I26[ a continuous increase of p as a function of /i signals a compressible 
phase, while a solid phase is characterized by a constant p for increasing fi. Panel 

(a) , corresponding to J/V = 0.05, shows a series of large constant-density plateaux 
connected by a progression of smaller steps and regions of continuous increase of p. 
Here, the main plateaux correspond to the Mott lobes of Fig. [25l while the other steps 
correspond to incompressible phases, with a fixed, integer, number of particles. This 
progression of steps is an indication of a Devil's-like staircase in the density, which 
was discussed in [13] for a one dimensional system. Instead, for J/V = 0.1 in panel 

(b) , no evidence of such a phase is found. 
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Appendix A. Parametrization 

Through the variation of (q, P), we aim at describing the transition between the states 
I initial and |$) final, as for example the one schematically represented in Fig. ITlT c). 
During this process there are sites initially occupied that become empty, like the blue- 
framed site of Fig. [TlT c). which we call the (B)-sites, and vice versa, like the site 
on the left of B, which is initially empty and occupied at the end, and we name the 
(A)-sites. When the initial and final states are non-degenerate, as for example the 
case sketched in Fig. [T2l d). there are also sites that do not change and remain either 
full (F) or empty (E). 

During this process, the Gutwiller amplitudes (|4.7[) have to be normalized at each 
site, and the total number of particles has to be conserved, namely 



where Ns is the total number of sites. We choose (g, P) = {qq, Pq) to be the variational 
parameters of the blue-framed site, and the normalization condition (|A.ip together 
with the conservation of the number of particles (jA.2p between A and B give us three 



(Al) 




(A2) 
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coupled equations 

q^-P^ + {qff - {P?f = 2 

{q^f-{P^f + {qtf-{Ptf = 2 (A.3) 

{<fif {prr + iitr in? - 2. 

As explained in [33], we make use of the following Ansatz 

qi =q 
= p 

\ (A.4) 

It is clear that a situation where site A is empty and site B is full corresponds to the 
values {q,P) = (0,0), while at {q, P) — (-^2,0) the contrary is true. For degenerate 
initial and final states as in the case of Fig. ITlT c). the remaining sites they either 
behave like A or B, which implies another set of conditions summarized as follows 

(i) A(B) 
% = % 

p(«) _ pA(B) 

" ' ' (A.5) 

p(i) „ pA(B) 
^1 — ^1 

depending on whether the site i is initially empty (A) or occupied (B). Instead, when 
the initial and final states are non-degenerate, as is the case considered in Fig. [T^d), 
there are also sites that we assume not to change and remain either full (F) or empty 
(E). For these sites the constraints are respectively given by 

ql -2 (A.6) 

Fo- = P! = 0, 

and 

qf =0 (A.7) 
P^ = Pf = 0. 

All these conditions (jA.3[) - (IA.7p . which we name Cc, enter explicitly into the calculation 
of Hamiltonian (|4.10p . 
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